Analyses script: EIH ~ Test modality and site
1 Library calls, set-up, data load
library(lattice) # dot plot for empty lme
library(moments) # distribution parameters
library(superb) # plot annotation
library(DescTools) # function Sign test for paired non-parametric test
library(effectsize) # effect sizes computations
library(grid) # plots organization
library(gridExtra) # plots organizaton
library(rstatix) # tidy base r stat
library(emmeans) # post-hoc tests
library(ggstatsplot) # plots extra functions
library(glue) # augmented paste function
library(tidyverse) # dataset organization tools
library(ftExtra) # package for tabular presentatino of results
library(officedown) # package for tabular presentatino of results
library(officer) # package for tabular presentatino of results
library(flextable) # package for tabular presentatino of results
library(rvg) # vectorial export of plots
library(patchwork) # plots organization
library(readxl) # import data from excel file
library(stringr) # function to work with strings
library(ggpubr) # added plotting functions
library(lme4) # linear mixed models computations
library(broom) # tidy model presentation
library(broom.mixed) # tidy model presentation for lme
library(ggpp) # added plotting functions
library(see) # added plotting functions
library(ggtext) # added plotting functions
library(ggsci) # added plotting functions
library(kableExtra) # html table
library(psych) # descriptive statistics and reliability assessment
library(here) # path definition
library(lmerTest) # compute p values for lme
library(ggplotify) # package allowing conversion of plot to ggplot objectsoptions(contrasts = c("contr.sum","contr.poly")) # setting a sum contrast/effect coding for lmer
set_flextable_defaults(font.size = 8,font.family = "Helevetica", padding=2, digits = 3) # display for flextables
knitr::opts_chunk$set( # display for chunks
tidy = TRUE,
warning = F,
comment = NULL,
message = F,
results = 'markup'
)
#fun_names <- dir(here("functions_wp1"))
##for(i in 1:length(fun_names)){
# source(here(paste0("functions_wp1/", fun_names[i])))
#}# confidence intervals
ci_low <- function(x){
avg <- mean(x, na.rm = T)
sd <- sd(x, na.rm = T)
n <- length(x)
error <- qnorm(0.95) * sd/sqrt(n)
avg - error
}
ci_high <- function(x){
avg <- mean(x, na.rm = T)
sd <- sd(x, na.rm = T)
n <- length(x)
error <- qnorm(0.95) * sd/sqrt(n)
avg + error
}
# mean function for reliab calculation
m1 <- function(x,y,z){
(x+y+z)/3
}
m2 <- function(x,y){
(x+y)/2
}
# mean function
m <- function(x,y,z){
(x+y+z)/3
}
## Parenthesis
par <- function(x,na.rm=FALSE){paste0("(",x,")")}
# formatting decimals
specify_decimal <- function(x, k) trimws(format(round(x, k), nsmall=k))
# formatting p values
pvalue_format <- function(x) {
z <- cut(x, breaks = c(-Inf, 0.001, 0.01, 0.05, 0.1, Inf), labels = c("***", " **", " *", " .", " "))
z <- as.character(z)
z[is.na(x)] <- ""
z
}
pvalue_format_table <- function(p){
if(p<0.001){
return("< .001")
}
else{
(return(format(p, digits = 3)))
}
}
pvalue_format_plot <- function(p){
if(p<0.001){
return("p < .001")
}
else{
(return(paste("p =",format(p, digits = 3))))
}
}
# ft formatting
dash_border <- officer::fp_border(color = "gray", style = "dashed")
# log t, will return NA if input is NA
logt <- function(x){
if(is.na(x)){
return(NA)
} else {
x <- abs(x) + 0.1
}
log <- log10(x)
return(log)
}
# function that returns the value at nth position from the max
# 1 is the max value
# input : vector, nth position
# output: value
tn <- function(x,n){
v <- sort(unique(x), decreasing = T)[n]
return(v)
}
# function to detect cooks values > 1 in a linear mixed model
cooks_outliers <- function(lmm, d, criteria = 1){
cooksd <- cooks.distance(lmm)
influential <- as.numeric(names(cooksd)[(cooksd > criteria)])
if(is_empty(influential)){
cat("No outlier detected")
} else {
cat(length(influential), "outliers detected. Observation(s) number: ", influential)
d.nouout <- d[-influential,]
return(d.nouout)
}
}
# ggplot custom theme
my_theme <- function(#base_size = 9, # theming function for ggplot
base_family = 'Roboto',
#base_line_size = base_size/22,
#base_rect_size = base_size/22,
border = TRUE) {
theme(axis.line.y = element_line(colour = "grey"),
axis.line.x = element_line(colour = "grey"),
panel.background = element_blank(),
plot.title = element_text(hjust = 0),
plot.subtitle = element_text(hjust = 0),
axis.title.y = element_markdown(),
plot.margin = unit(c(0, 0.5, 0.5, 0), "cm"))
}
var_summary <- function(data, var_vec){
# summary of non-grouped data
smry <- data %>%
select({{var_vec}}) %>%
psych::describe(quant=c(.25,.75), IQR = T) %>%
as_tibble(rownames="rowname") %>%
select(var=rowname, mean, sd, min, median, max, IQR) |>
mutate(across(is.numeric,round, 2))
return(smry)
}
# Vladimir Aron 22/07/2024
# input: data set, title, units, df anova, var2 = variable to filter anova table with
# output: an interaction plot with the results of the interaction computed from the anova df using corr/non-corr p value
plot_prepost_an <- function(d, t, u, df_an,v, y = "test_value", x = T, corr = F){
# d <- dm.ppt.quad
# t <- "CDT"
# u <- "°C"
# mod <- m.0.pptquad
# y <- "test_value"
# df_an <- an_mpp
# v = "CDT_Quadriceps"
# d = dm.adt
# t = "Auditory detection threshold"
# u = "mA"
# df_an = an_mpp
# v = "ADT_No site"
main <- strsplit(v, "_")[[1]][1]
sub <- strsplit(v, "_")[[1]][2]
d <- d |>
mutate(timing = factor(timing, levels = c("pre", "post"), labels = c("Pre", "Post")),
condition = factor(condition, levels = c("cont", "exp"), labels = c("Control", "Exercise")))
d.sum <- d |>
group_by(timing, condition, ID) |>
mutate(tm = mean(.data[[y]], na.rm = T)) |>
slice(1) |>
ungroup() |>
group_by(timing, condition) |>
summarise(m = mean(tm, na.rm = T),
cil = ci_low(tm),
cih = ci_high(tm))
d.all <- d |>
group_by(timing, condition, ID) |>
mutate(tm = mean(.data[[y]], na.rm = T)) |>
slice(1) |>
mutate(ID2 = factor(paste(condition, ID))) |>
arrange(ID2)
grand.mean <- mean(d.all$tm, na.rm = T)
if(corr == T){
label <- df_an |>
filter(var2 == v) |>
filter(term == "timing:condition") |>
mutate(p.corr = pvalue_format_plot(p.corr)) |>
mutate(interaction = glue("Timing x Condition: {p.corr}")) |>
pull(interaction)
} else{
label <- df_an |>
filter(var2 == v) |>
filter(term == "timing:condition") |>
mutate(p.value = pvalue_format_plot(p.value)) |>
mutate(interaction = glue("Timing x Condition: {p.value}")) |>
pull(interaction)
}
grob <- grobTree(textGrob(label, x = 0.1, y = 0.95, hjust=0,
gp=gpar(col="black", fontsize=9, fontfamily = "Roboto")))
grob2 <- grobTree(textGrob(label, x = 0.1, y = 0.95, hjust=0,
gp=gpar(col="black", fontsize=9, fontfamily = "Roboto")))
pos <- position_dodge2(width = 0.1)
pos2 <- position_dodgenudge(width = 0.3)
# pos <- position_dodgenudge(width = 0.3)
pos3 <- position_dodge2(width = 0.2)
p_full <- d.sum |>
ggplot(aes(x = timing, y = m, col = condition)) +
geom_point(position = pos2, size = 3) +
geom_line(aes(group = factor(condition)), position = pos2, linewidth = 1) +
geom_errorbar(aes(ymin =cil, ymax = cih), width = 0.1, position = pos2) +
geom_point(data = d.all, aes(y = tm, x = timing, col = condition) ,alpha = 0.2, size = 1, position = pos) +
geom_line(data = d.all, aes(y = tm, x = timing, col = condition, group = ID2) ,alpha = 0.05, position = pos) +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
my_theme() +
labs(title = t,
x = "Timing",
col = "Condition",
y = u) + annotation_custom(grob)
if(main == "ADT"){
coord_y_lim = c(quantile(d.all$tm, 0.1, na.rm = T), quantile(d.all$tm, 0.99, na.rm = T))
} else {
coord_y_lim = c(quantile(d.all$tm, 0.2), quantile(d.all$tm, 0.8))
}
p_means <- d.sum |>
ggplot(aes(x = timing, y = m, col = condition)) +
geom_point(position = pos2, size = 4, show.legend = F) +
geom_line(aes(group = factor(condition)), position = pos2, linewidth = 1, show.legend = F) +
geom_errorbar(aes(ymin =cil, ymax = cih), width = 0.05, position = pos2, show.legend = F) +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(aes(linetype = "Overall mean", yintercept = grand.mean), color = "lightgrey") +
scale_linetype_manual(values = c("Overall mean" = "dashed")) +
coord_cartesian(ylim = coord_y_lim) +
my_theme() +
labs(
x = "Timing",
col = "Condition",
linetype = "",
y = paste("Condition means (", u, ")")) + annotation_custom(grob2)+
theme(legend.position = "top")
p_points <- ggplot()+
geom_point(data = d.all, aes(y = tm, x = timing, col = condition) , position = pos3, alpha = 0.7, size = 2, show.legend = F) +
geom_line(data = d.all, aes(y = tm, x = timing, col = condition, group = ID2) ,alpha = 0.2, position = pos3, linewidth = 0.5, show.legend = F) +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
my_theme() +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
theme(legend.position = "none") +
geom_hline(aes(linetype = "Overall mean", yintercept = grand.mean), color = "lightgrey", show.legend = F) +
scale_linetype_manual(values = c("Overall mean" = "dashed")) +
labs(title = t,
subtitle = sub,
x = "Timing",
linetype = "",
y = paste("Individual means (", u, ")")) +
theme(legend.position = "none")
grob_cont <- grobTree(textGrob("Control", x = 0.05, y = 0.95, hjust=0,
gp=gpar(col="#00468BFF", fontsize=9, fontfamily = "Roboto")))
grob_exp <- grobTree(textGrob("Exercise", x = 0.05, y = 0.95, hjust=0,
gp=gpar(col="#ED0000FF", fontsize=9, fontfamily = "Roboto")))
pos4 <- position_dodge2(width = 0.1)
p_points_cont <- d.all |>
filter(condition == "Control") |>
ggplot()+
geom_point(aes(y = tm, x = timing, col = condition) , position = pos4, alpha = 0.7, size = 2, show.legend = F) +
geom_line(aes(y = tm, x = timing, col = condition, group = ID2) ,alpha = 0.2, linewidth = 0.5, position = pos4, show.legend = F) +
scale_color_manual(values = "#00468BFF") +
my_theme() +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(aes(linetype = "Overall mean", yintercept = grand.mean), color = "lightgrey", show.legend = F) +
scale_linetype_manual(values = c("Overall mean" = "dashed")) +
labs(title = t,
subtitle = sub,
x = "Timing",
y = paste("Individual means (", u, ")") ) +
theme(
legend.position = "none"
) + annotation_custom(grob_cont)
p_points_exp <- d.all |>
filter(condition == "Exercise") |>
ggplot()+
geom_point(aes(y = tm, x = timing, col = condition) , position = pos4, alpha = 0.7, size = 2, show.legend = F) +
geom_line(aes(y = tm, x = timing, col = condition, group = ID2) ,alpha = 0.2, linewidth = 0.5, position = pos4, show.legend = F) +
scale_color_manual(values = "#ED0000FF") +
my_theme() +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(aes(linetype = "Overall mean", yintercept = grand.mean), color = "lightgrey", show.legend = F) +
scale_linetype_manual(values = c("Overall mean" = "dashed")) +
labs(title = "",
subtitle = "",
x = "Timing",
y ="" ) + annotation_custom(grob_exp)
if(x == F){
p_points_cont <- p_points_cont + theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()) +
labs(x = "")
p_points_exp <- p_points_exp + theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()) +
labs(x = "")
p_points <- p_points +
scale_x_discrete(breaks = "") +
labs(x = "")
p_means <- p_means +
scale_x_discrete(breaks = "") +
labs(x = "")
}
# p <- p_points + p_means
p <- (p_points_cont + p_points_exp + p_means)
#annotate("text", x = 1.2, y = max(d.all$tm), label = label)
if(y == "test_value"){
return(p)
} else {
return(p_full)
}
}
# Vladimir Aron 22/07/2024
# input: data set, title, units, lmm model
# output : interaction plot site:condition + annotation of the significance of the interaction
plot_eih_s <- function(d,t, u, mod, y_axis = T, x = F){
#d <- dm.ppt
# mod <- m.0.ppt
# t <- "d"
# u <- "d"
label <- anova(mod, ddf = "Kenward-Roger") |>
tidy() |>
filter(term == "condition:site") |>
mutate(p = pvalue_format_plot(p.value)) |>
mutate(interaction = glue("Site x Condition: {p}")) |>
pull(interaction)
grob <- grobTree(textGrob(label, x = 0.05, y = 0.95, hjust=0,
gp=gpar(col="black", fontsize=9, fontfamily = "Roboto")))
grob_cont <- grobTree(textGrob("Control", x = 0.05, y = 0.95, hjust=0,
gp=gpar(col="#00468BFF", fontsize=9, fontfamily = "Roboto")))
grob_exp <- grobTree(textGrob("Exercise", x = 0.05, y = 0.95, hjust=0,
gp=gpar(col="#ED0000FF", fontsize=9, fontfamily = "Roboto")))
dd <- d |>
mutate(site = factor(site, levels = c("forearm", "quad"), labels = c("Forearm", "Quadriceps")),
condition = factor(condition, levels = c("cont", "exp"), labels = c("Control", "Exercise")))
d.sum <- dd |>
group_by(site, condition, ID) |>
mutate(tm = mean(abc, na.rm = T)) |>
slice(1) |>
ungroup() |>
group_by(site, condition) |>
summarise(m = mean(tm),
cil = ci_low(tm),
cih = ci_high(tm))
d.all <- dd |>
group_by(site, condition, ID) |>
mutate(tm = mean(abc, na.rm = T)) |>
slice(1) |>
mutate(ID2 = factor(paste(ID, condition))) |>
arrange(ID2)
grand.mean <- mean(d.all$tm, na.rm = T)
coord_y_lim = c(quantile(d.all$tm, 0.1), quantile(d.all$tm, 0.9))
pos <- position_dodge2(width = 0.1)
pos2 <- position_dodgenudge(width = 0.3)
p.all <- d.sum |>
ggplot(aes(x = site, y = m, col = condition, group = condition)) +
geom_point(position = pos2, size = 4) +
geom_line(position = pos2, linewidth = 1) +
geom_errorbar(aes(ymin =cil, ymax = cih), width = 0.1, position = pos2) +
geom_point(data = d.all, aes(y = tm, x = site, col = condition) ,alpha = 0.2, size = 1, position = pos) +
geom_line(data = d.all, aes(y = tm, x = site, col = condition, group = ID2) ,alpha = 0.05, position = pos) +
scale_color_manual(values = c("#ED0000FF","#00468BFF")) +
geom_hline(yintercept = 0, col = "grey", alpha = 0.5, linetype = "dashed") +
my_theme() +
labs(title = t,
x = "Site",
col = "Condition",
y = u) +
theme(legend.position = "top")+ annotation_custom(grob)
if(y_axis == T){
# y_title_points <- paste("Diff<sub>post-pre</sub> (", u, ")")
# y_title_means <- paste("Diff<sub>post-pre</sub> (", u, ")")
y_title_points <- y_title_means <- "Diff<sub>post-pre</sub> (standardized units)"
}else{
y_title_points <- y_title_means <- ""
}
pos3 <- position_dodge2(width = 0.2)
p_points <- ggplot()+
geom_point(data = d.all, aes(y = tm, x = site, col = condition) , position = pos3, alpha = 0.7, size = 2, show.legend = F) +
geom_line(data = d.all, aes(y = tm, x = site, col = condition, group = ID2) ,alpha = 0.2, linewidth = 0.5, position = pos3, show.legend = F) +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
my_theme() +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(yintercept = 0, col = "grey", alpha = 0.5, linetype = "dashed") +
labs(title = t,
subtitle = "Individual means",
x = "Site",
y = y_title_points ) +
theme(
legend.position = "none"
)
pos4 <- position_dodge2(width = 0.1)
p_points_cont <- d.all |>
filter(condition == "Control") |>
ggplot()+
geom_point(aes(y = tm, x = site, col = condition) , position = pos4, alpha = 0.7, size = 2, show.legend = F) +
geom_line(aes(y = tm, x = site, col = condition, group = ID2) ,alpha = 0.2, linewidth = 0.5, position = pos4, show.legend = F) +
scale_color_manual(values = "#00468BFF") +
my_theme() +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
labs(title = t,
subtitle = "Individual means",
x = "Site",
y = y_title_points ) +
theme(
legend.position = "none"
) + annotation_custom(grob_cont)
p_points_exp <- d.all |>
filter(condition == "Exercise") |>
ggplot()+
geom_point(aes(y = tm, x = site, col = condition) , position = pos4, alpha = 0.7, size = 2, show.legend = F) +
geom_line(aes(y = tm, x = site, col = condition, group = ID2) ,alpha = 0.2, linewidth = 0.5, position = pos4, show.legend = F) +
scale_color_manual(values = "#ED0000FF") +
my_theme() +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
labs(title = "",
subtitle = "",
x = "Site",
y ="" ) +
theme(
legend.position = "none"
) + annotation_custom(grob_exp)
p_means <- d.sum |>
ggplot(aes(x = site, y = m, col = condition)) +
geom_point(position = pos2, size = 4) +
geom_line(aes(group = factor(condition)), position = pos2, linewidth = 1) +
geom_errorbar(aes(ymin =cil, ymax = cih), width = 0.05, position = pos2) +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
scale_y_continuous(limits = c(min(d.all$tm), max(d.all$tm))) +
geom_hline(yintercept = 0, col = "black", linetype = "dashed") +
coord_cartesian(ylim = coord_y_lim) +
my_theme() +
labs(
subtitle = "Condition means",
col = "Condition",
x = "Site",
y = y_title_means) +
annotation_custom(grob)
if(x == F){
p_points_cont <- p_points_cont + theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()) +
labs(x = "")
p_points_exp <- p_points_exp + theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()) +
labs(x = "")
p_points <- p_points + theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()) +
labs(x = "")
p_means <- p_means + theme(axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()) +
labs(x = "")
} else{
p_points <- p_points
p_means <- p_means
}
p_points_sep <- p_points_cont + p_points_exp
res <- list(
p.all = p.all,
p_ponts_sep = p_points_sep,
p_means = p_means,
p_points = p_points
)
#p <- p +
# annotate("text", x = 1, y = max(d.all$tm), label = label)
return(res)
}# data
n <- 40 # enter number of desired participants to analyse
excluded <- c("35", "40")
n_an <- n - length(excluded)
# data screening session
## Data baseline
data_b <- read_excel(here("crf_eih.xlsx"), na = c("NM","missing"), sheet = "baseline", skip = 1, col_names = TRUE) %>% select(!dominance)
data_b <- data_b[1:n,]
## data test session
data_t <- read_excel(here("crf_eih.xlsx"), na = c("NM","missing"), sheet = "test_data", skip = 0, col_names = TRUE)
data_t <- data_t[1:n,]
## merging data sets
data_bt <- full_join(data_b, data_t)
# data sessions
data_1 <- read_excel(here("crf_eih.xlsx"), na = c("NM","missing"), sheet = "session_1", skip = 2, col_names = TRUE) # import data from session 1
data_1 <- data_1[1:n,]
data_2 <- read_excel(here("crf_eih.xlsx"), na = c("NM","missing"), sheet = "session_2", skip = 2, col_names = TRUE) # import data from session 2
data_2 <- data_2[1:n,]
data_eih <- rbind(data_1, data_2) # combine data from session 1 and 2
# merging all data sets
data_tot <- full_join(data_eih, data_bt)
data <- data_tot %>%
filter(!ID %in% excluded) |>
mutate(IPAQ_category = factor(IPAQ_category, levels = c("LOW", "MODERATE", "HIGH"))) |>
rename(c("SSS" = "ESS", "RPE" = "rpe_post"))
# removing values ADT post for subject 10 --> technical issue (script ran but somehow reloading it reduced the values to more "normal" levels )
data["adt_post_1"][data["adt_post_1"] == 0.46] <- "NA"
data["adt_post_2"][data["adt_post_2"] == 0.41] <- "NA"
data["adt_post_3"][data["adt_post_3"] == 0.43] <- "NA"
#data |>
# filter(condition == "exp") |>
#select(ID, temp_forearm_post1_nd , temp_forearm_post2_nd) |> view()
#save(data, "data_testmod_site_20240814")# data sets are put in long format
d_tms_long <- data %>% # data of test mod for each site, condition, order and timing put in long format
mutate(across(starts_with("ppt_"),as.numeric)) %>%
mutate(across(starts_with("cdt_"),as.numeric)) %>%
mutate(across(starts_with("hpt_"),as.numeric)) %>%
mutate(across(starts_with("pp_"),as.numeric)) %>%
mutate(across(starts_with("adt_"),as.numeric)) %>%
rename_with(~sub("adt_", "adt_nosite_", .), everything()) %>%
select(ID, condition, session, starts_with(c("ppt_", "cdt_", "hpt_", "pp_", "adt_"))) |>
pivot_longer(cols = starts_with(c("ppt_", "cdt_", "hpt_", "pp_", "adt_")),
names_to = c("test_mod","site", "timing", "order"),
names_pattern = "(^\\w+)_(\\w+)_(\\w+)_(\\w+)",
values_to = c("test_value")
)
reference_date <- as.POSIXct("1970-01-01 00:00:00", tz = "UTC")
d_cov_long <- data %>% # data set of covariables (blood pressure, heart rate...) is put in long format
mutate(across(starts_with(c("temp_q", "temp_f")),as.numeric)) %>%
mutate(across(starts_with("hr_"),as.numeric)) %>%
mutate(across(starts_with("br_"),as.numeric)) %>%
rowwise() |>
# bp
separate(BP_rest, c("sbp_pre", "dbp_pre"), sep = "/") %>%
separate(bp_post, c("sbp_post", "dbp_post"), sep = "/") %>%
mutate(across(c(starts_with("sbp_") , starts_with("dbp_")), as.numeric)) |>
mutate(slope_bp_sys = sbp_post - sbp_pre,
slope_bp_diast = dbp_post - dbp_pre) %>% # absolute change of blood pressure
select(ID, condition,HR_rest, BR_rest, duration_pre, hr_post, sbp_post, dbp_post, br_post, duration_post, sbp_pre, dbp_pre) |>
rename(c(hr_pre = HR_rest, br_pre = BR_rest)) |>
group_by(ID, condition) |>
slice(1) |>
mutate(duration_pre = str_replace(duration_pre, ";", ":"),
duration_post = str_replace(duration_post, ";", ":")) %>%
mutate(duration_pre = strptime(duration_pre, "%M:%S"),
duration_post = strptime(duration_post, "%M:%S")) %>%
mutate(duration_pre = seconds_to_period(duration_pre),
duration_post = seconds_to_period(duration_post)) %>%
mutate(duration_pre = period_to_seconds(duration_pre),
duration_post = period_to_seconds(duration_post)) |>
pivot_longer(cols = c(everything(), -ID,-condition ), names_to = c("var","timing"),names_sep = "_", values_to = "values" ) |>
pivot_wider(names_from = var, values_from = values) |>
mutate(duration = as.numeric(seconds_to_period(duration)) + reference_date) |>
mutate(duration = minute(duration) + second(duration)/60)
d_mean_temp_long <- data |> # skin temperature data are averaged before and after conditions and put in long format
select(ID,condition, starts_with("temp")) |>
select(-temp_room) |>
mutate(across(starts_with("temp"), as.numeric)) |>
rowwise() |>
mutate(temp_quad_pre = (temp_quad_pre1_d + temp_quad_pre2_d)/2,
temp_forearm_pre = (temp_forearm_pre1_nd + temp_forearm_pre2_nd)/2,
temp_quad_post = (temp_quad_post1_d + temp_quad_post2_d)/2,
temp_forearm_post = (temp_forearm_post1_nd + temp_forearm_post2_nd)/2,
) |>
select(ID, condition, temp_quad_pre, temp_forearm_pre, temp_quad_post, temp_forearm_post) |>
group_by(ID) |>
pivot_longer(cols = starts_with("temp"),
names_to = c("xx", "site", "timing"),
names_sep = "_",
values_to = "temp") |>
select(-xx)
data_long <- full_join(d_tms_long, d_mean_temp_long) |> # data set for skin temperature and tms are merged
mutate(site = if_else(test_mod == "adt", "no site", site)) |>
mutate(session = factor(session, levels = c("session_1", "session_2")),
timing = factor(timing, levels = c("pre", "post")),
order = factor(order, levels = c("1", "2", "3")),
test_mod = factor(test_mod, levels = c("ppt", "pp", "hpt", "cdt", "adt")),
site = factor(site, levels = c("forearm", "quad", "no site")),
condition = factor(condition, levels = c("exp", "cont"))
) |>
mutate(test_value = if_else(test_mod == "cdt",abs(test_value - 32), test_value)) |> # conversion of CDT in absolute difference from baseline
mutate(test_value = if_else(test_mod == "pp", 105 - test_value, test_value)) |> # reversing scale of PP
mutate(test_value = if_else(test_mod == "pp",(test_value/105)*100 ,test_value ))|> # data were collected on a 105mm scale > they are converted to a 100mm scale
mutate(test_value = if_else(test_mod == "ppt", test_value*10, test_value)) # converting N/cm2 in kilospascale
check_mod <-function(lmm,data){
mod_name <- deparse(substitute(lmm))
data_name <- deparse(substitute(data))
mod.info <- paste(mod_name, data_name, sep = ", ")
data[,"predicted"] <- predict(lmm)
data[,"resid"] <- residuals(lmm)
#fitted residuals
pfit <- ggplot(data = {{data}}, aes(x = predicted, y = resid)) +
geom_point() +
geom_hline(yintercept = 0, colour = "red") +
facet_wrap(site~test_mod , scales = "free") +
labs(title = paste("Fitted values ~ Residuals; ", mod.info)) + my_theme()
#plot(fitted(lmm),residuals(lmm), main = paste("Fitted values ~ Residuals; ", mod.info))
#normalité epsilon
pqq <- ggplot(data = data, aes(sample = resid)) +
stat_qq() +
stat_qq_line(color = "red") +
facet_wrap(site~test_mod, scales = "free") +
labs(title = paste("Normality residuals; ", mod.info)) + my_theme()
phist <- ggplot(data = data, aes(x = resid)) +
geom_histogram(aes(y = after_stat(density)), bins = 9) +
geom_density(color = "red") +
facet_wrap(site~test_mod, scales = "free") + my_theme() +
labs(title = paste("Normality residuals; ", mod.info))
# qqnorm(residuals(lmm), main = paste("Normality residuals; ", mod.info) )
#qqline(residuals(lmm))
# qqplot normalite u_0j
raneff <- names(ranef(lmm))
if (length(raneff) > 1){
cat("There are more than one random effets: ", paste(raneff, collapse = ", "), "\n")
ran_name <- readline(prompt = "Enter name of desired random effet: ")
ran_num <- which(raneff == ran_name)
str(ranef(lmm))
ranef(lmm)[[ran_num]][,1]
df <- data.frame(ran = ranef(lmm))
qqnorm(ranef(lmm)[[ran_num]][,1], main = paste("Normality random effects: ", ran_name, ", ", mod.info))
qqline(ranef(lmm)[[ran_num]][,1], main = paste("Normality random effects: ", ran_name, ", ", mod.info))
} else {
qqnorm(ranef(lmm)$ID[,1], main = paste("Normality random effects; ", mod.info))
qqline(ranef(lmm)$ID[,1])
}
pcook <- ggplot(data.frame(cook=cooks.distance(lmm),id=data$ID),
aes(x=cook,y=id)) +
geom_point() +
labs(title = paste("ID ~ Cook's distance; ", mod.info)) +
theme_bw()
#(pfit + pcook)/(pqq + phist)
print(pfit)
print(pcook)
print(pqq)
print(phist)
}
check_mod_tms <-function(lmm,data){
mod_name <- deparse(substitute(lmm))
data_name <- deparse(substitute(data))
mod.info <- paste(mod_name, data_name, sep = ", ")
data[,"predicted"] <- predict(lmm)
data[,"resid"] <- residuals(lmm)
#fitted residuals
pfit <- ggplot(data = {{data}}, aes(x = predicted, y = resid)) +
geom_point() +
geom_hline(yintercept = 0, colour = "red") +
facet_wrap(timing~condition , scales = "free") +
labs(title = paste("Fitted values ~ Residuals; ", mod.info)) + my_theme()
#plot(fitted(lmm),residuals(lmm), main = paste("Fitted values ~ Residuals; ", mod.info))
#normalité epsilon
pqq <- ggplot(data = data, aes(sample = resid)) +
stat_qq() +
stat_qq_line(color = "red") +
facet_wrap(timing~condition, scales = "free") +
labs(title = paste("Normality residuals; ", mod.info)) + my_theme()
phist <- ggplot(data = data, aes(x = resid)) +
geom_histogram(aes(y = after_stat(density)), bins = 9) +
geom_density(color = "red") +
facet_wrap(timing~condition, scales = "free") + my_theme() +
labs(title = paste("Normality residuals; ", mod.info))
# qqnorm(residuals(lmm), main = paste("Normality residuals; ", mod.info) )
#qqline(residuals(lmm))
# qqplot normalite u_0j
raneff <- names(ranef(lmm))
if (length(raneff) > 1){
cat("There are more than one random effets: ", paste(raneff, collapse = ", "), "\n")
ran_name <- readline(prompt = "Enter name of desired random effet: ")
ran_num <- which(raneff == ran_name)
str(ranef(lmm))
ranef(lmm)[[ran_num]][,1]
df <- data.frame(ran = ranef(lmm))
qqnorm(ranef(lmm)[[ran_num]][,1], main = paste("Normality random effects: ", ran_name, ", ", mod.info))
qqline(ranef(lmm)[[ran_num]][,1], main = paste("Normality random effects: ", ran_name, ", ", mod.info))
} else {
qqnorm(ranef(lmm)$ID[,1], main = paste("Normality random effects; ", mod.info))
qqline(ranef(lmm)$ID[,1])
}
pcook <- ggplot(data.frame(cook=cooks.distance(lmm),id=data$ID),
aes(x=cook,y=id)) +
geom_point() +
labs(title = paste("ID ~ Cook's distance; ", mod.info)) +
theme_bw()
#(pfit + pcook)/(pqq + phist)
print(pfit)
print(pcook)
print(pqq)
print(phist)
}
check_mod_s <-function(lmm,data){
mod_name <- deparse(substitute(lmm))
data_name <- deparse(substitute(data))
mod.info <- paste(mod_name, data_name, sep = ", ")
data[,"predicted"] <- predict(lmm)
data[,"resid"] <- residuals(lmm)
#fitted residuals
pfit <- ggplot(data = {{data}}, aes(x = predicted, y = resid)) +
geom_point() +
geom_hline(yintercept = 0, colour = "red") +
facet_wrap(site~condition , scales = "free") +
labs(title = paste("Fitted values ~ Residuals; ", mod.info)) + my_theme()
#plot(fitted(lmm),residuals(lmm), main = paste("Fitted values ~ Residuals; ", mod.info))
#normalité epsilon
pqq <- ggplot(data = data, aes(sample = resid)) +
stat_qq() +
stat_qq_line(color = "red") +
facet_wrap(site~condition, scales = "free") +
labs(title = paste("Normality residuals; ", mod.info)) + my_theme()
phist <- ggplot(data = data, aes(x = resid)) +
geom_histogram(aes(y = after_stat(density)), bins = 9) +
geom_density(color = "red") +
facet_wrap(site~condition, scales = "free") + my_theme() +
labs(title = paste("Normality residuals; ", mod.info))
# qqnorm(residuals(lmm), main = paste("Normality residuals; ", mod.info) )
#qqline(residuals(lmm))
# qqplot normalite u_0j
raneff <- names(ranef(lmm))
if (length(raneff) > 1){
cat("There are more than one random effets: ", paste(raneff, collapse = ", "), "\n")
ran_name <- readline(prompt = "Enter name of desired random effet: ")
ran_num <- which(raneff == ran_name)
str(ranef(lmm))
ranef(lmm)[[ran_num]][,1]
df <- data.frame(ran = ranef(lmm))
qqnorm(ranef(lmm)[[ran_num]][,1], main = paste("Normality random effects: ", ran_name, ", ", mod.info))
qqline(ranef(lmm)[[ran_num]][,1], main = paste("Normality random effects: ", ran_name, ", ", mod.info))
} else {
qqnorm(ranef(lmm)$ID[,1], main = paste("Normality random effects; ", mod.info))
qqline(ranef(lmm)$ID[,1])
}
pcook <- ggplot(data.frame(cook=cooks.distance(lmm),id=data$ID),
aes(x=cook,y=id)) +
geom_point() +
labs(title = paste("ID ~ Cook's distance; ", mod.info)) +
theme_bw()
#(pfit + pcook)/(pqq + phist)
print(pfit)
print(pcook)
print(pqq)
print(phist)
}2 Participants data
Summary statistics of participants data (age, bmi, laterality, and level of education) as well as scores to questionnaires are computed and formatted in a table.
var_new_names <- c("age" = "Age (years)" ,"bmi" = "BMI (kg/m^2^)" ,"FlANDERS_lat"= "FLANDERS (n)" , "level_education" = "Level of education (n)" , "W_peak_rel" = "Relative PPO (W/kg)", "HR_max" = "Maximal HR (bpm)" , "STAI_2"= "STAI-Y2 (/80)" , "PCS" = "PCS total (/51)" , "FPQ" = "FPQ-15 (/100)" , "IPAQ_category" = "IPAQ category (n)", "IPAQ_total_moderate"= "IPAQ total mPA (MET/min/wk)" , "IPAQ_total_vigourous" = "IPAQ total vPA (MET/min/wk)" , "IPAQ_total_pa"= "IPAQ total PA (MET/min/wk)", "W_peak" = "PPO (W)", "TSK" = "TSK (/37)" )
cont_data <- var_summary(data, var_vec = c(age, bmi, W_peak, W_peak_rel, HR_max, STAI_2, PCS, FPQ, TSK, IPAQ_total_moderate, IPAQ_total_vigourous, IPAQ_total_pa)) |>
mutate(value_1 = glue("{mean} ± {sd}")) |>
mutate(value_2 = glue("{median} ({IQR})")) |>
mutate(value = if_else(var %in% c("IPAQ_total_moderate", "IPAQ_total_vigourous", "IPAQ_total_pa"), value_2, value_1)) |>
select(var, value)
cont_data$levels <- c(rep(NA, nrow(cont_data)))
ipaq_cat <- as_tibble(table(subset(data, condition == "exp", IPAQ_category))) |>
rename(levels = IPAQ_category) |>
rename(value = n)
ipaq_cat$var <- c(rep("IPAQ_category",nrow(ipaq_cat)))
flanders <- as_tibble(table(subset(data, condition == "exp", FlANDERS_lat))) |>
rename(levels = FlANDERS_lat) |>
rename(value = n)
flanders$var <- c(rep("FlANDERS_lat",nrow(flanders)))
educ <- as_tibble(table(subset(data, condition == "exp", level_education))) |>
rename(levels = level_education) |>
rename(value = n)
educ$var <- c(rep("level_education",nrow(educ)))
order_rows <- c("age", "bmi", "FLANDERS_lat", "level_education", "W_peak", "W_peak_rel", "HR_max", "STAI_2", "PCS", "FPQ", "TSK", "IPAQ_total_moderate", "IPAQ_total_vigourous", "IPAQ_total_pa", "IPAQ_category")
t1 <- rbind(cont_data, flanders, ipaq_cat) |>
select(var, levels, value) |>
arrange(factor(var, levels = order_rows)) |>
mutate(var = str_replace_all(var, var_new_names)) |>
flextable() %>%
set_header_labels(var = "Variables",levels = " ", value = paste("n =", n_an)) %>%
bold(part ="header") %>%
merge_v(j = 1, part = "body") |>
autofit() %>%
fix_border_issues() %>%
colformat_md(part = "body", j = 1) |>
hline(i = c(1:12,15), border = dash_border) |>
add_footer_lines(values = c("Data presented as mean ± SD; except for as median (interquartile range). \n Abbreviations: BMI, Body mass Index; PPO, Peak power output; HR, Heart rate; STAI-Y2, State-trait anxiety inventory (form 2, trait); PCS, Pain catastrophizing scale; FPQ-15; Fear of pain questionnaire 15 items; TSK, Tampa scale for kinesiophobia; IPAQ, International physical activity questionnaire; FLANDERS, Flinders handedness survey. "
)) %>%
#align(i = 1, align = "left", part = "footer") %>%
valign(part = "body", valign = "top") |>
set_table_properties(layout = "autofit")
#print(ts1, preview = "docx")
t1Variables |
| n = 38 |
|---|---|---|
Age (years) | 23.84 ± 3.69 | |
BMI (kg/m2) | 22.37 ± 2.63 | |
PPO (W) | 258.11 ± 44.27 | |
Relative PPO (W/kg) | 3.61 ± 0.64 | |
Maximal HR (bpm) | 184.68 ± 9.33 | |
STAI-Y2 (/80) | 34.32 ± 7.46 | |
PCS total (/51) | 13.32 ± 7.12 | |
FPQ-15 (/100) | 35.74 ± 12.08 | |
TSK (/37) | 32.82 ± 6.29 | |
IPAQ total mPA (MET/min/wk) | 1080 (1530) | |
IPAQ total vPA (MET/min/wk) | 2880 (2880) | |
IPAQ total PA (MET/min/wk) | 6339 (4974) | |
IPAQ category (n) | LOW | 1 |
MODERATE | 5 | |
HIGH | 32 | |
FLANDERS (n) | Ambidextrous | 1 |
Left | 3 | |
Right | 34 | |
Data presented as mean ± SD; except for as median (interquartile range). | ||
save_as_docx(t1, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/particpants_data.docx")3 Manipulation check
3.1 Control vs Exercise condition
Control vs Exercise conditions are compared, before and after intervention, on heart rate (HR), breathing rate (BR), and blood pressure (systolic, SBP; diastolic, DBP) . Scores of rate of perceived exertion (RPE), stanford sleepiness scale (SSS), state anxiety (STAI-Y1) and room temperature and humidity are also compared.
# function that returns a data set for a covariable
# data is groupd by ID, timing, condition
#input, var
# output, data set
d.mod.cov <- function(var){
dm <- d_cov_long |>
group_by(ID, timing, condition) |>
slice(1) |>
select({{var}}, timing, condition, ID) |>
rename(value = {{var}}) |>
mutate(timing = factor(timing, levels = c("pre", "post")),
condition = factor(condition, levels = c("cont", "exp") ))
dsum <- dm |>
ungroup() |>
group_by(timing, condition) |>
summarise(m = mean(value, na.rm = T),
cil = ci_low(value),
cih = ci_high(value))
res <- list(
dm = dm,
dsum = dsum
)
return(res)
}3.1.1 HR
d.hr <- d.mod.cov(hr)$dm
mod.hr <- lmer(value ~ timing*condition + (1|ID), data = d.hr)
cooks_outliers(mod.hr, d.hr)No outlier detectedcheck_mod_tms(mod.hr, d.hr)
The random effects follow a normal distribution with fat tails, relatively symmetrical. Residuals are homoskedastic. Residuals do not follow normality in the post:exp factor. Cooks distance < 1 suggest little influence of deviated data point.
an_hr <- anova(mod.hr, ddf = "Kenward-Roger") |> tidy()
an_hr |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| timing | 93307.60 | 93307.60 | 1 | 111 | 2179.184 | 0 |
| condition | 68467.60 | 68467.60 | 1 | 111 | 1599.050 | 0 |
| timing:condition | 68722.53 | 68722.53 | 1 | 111 | 1605.003 | 0 |
3.1.2 BR
# BR
d.br <- d.mod.cov(br)$dm
mod.br <- lmer(value ~ timing*condition + (1|ID), data = d.br)
cooks_outliers(mod.br, d.br)No outlier detectedcheck_mod_tms(mod.br, d.br)
Normal effects appear slighlty deviated from normality. Residuals are homoskedastic. Residuals do not follow normality in the post:exp factor. Cooks distance < 1 suggest little influence of deviated data point.
an_br <- anova(mod.br, ddf = "Kenward-Roger") |> tidy()
an_br |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| timing | 5887.605 | 5887.605 | 1 | 111 | 138.522 | 0 |
| condition | 3316.447 | 3316.447 | 1 | 111 | 78.028 | 0 |
| timing:condition | 3467.605 | 3467.605 | 1 | 111 | 81.585 | 0 |
3.1.3 SBP
d.sbp <- d.mod.cov(sbp)$dm
mod.sbp <- lmer(value ~ timing*condition + (1|ID), data = d.sbp)
cooks_outliers(mod.sbp, d.sbp)No outlier detectedcheck_mod_tms(mod.sbp, d.sbp)
Random effects are normally distributed. Residuals are homoskedastic and normally distributed across factors. No outlier is detected.
an_sbp <- anova(mod.sbp, ddf = "Kenward-Roger") |> tidy()
an_sbp |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| timing | 3720.421 | 3720.421 | 1 | 111 | 98.744 | 0 |
| condition | 3335.158 | 3335.158 | 1 | 111 | 88.519 | 0 |
| timing:condition | 2480.237 | 2480.237 | 1 | 111 | 65.828 | 0 |
3.1.4 DBP
d.dbp <- d.mod.cov(dbp)$dm
mod.dbp <- lmer(value ~ timing*condition + (1|ID), data = d.dbp)
d.dbp.noOut <- cooks_outliers(mod.dbp, d.dbp)1 outliers detected. Observation(s) number: 98mod.dbp.noOut <- lmer(value ~ timing*condition + (1|ID), data = d.dbp.noOut)
check_mod_tms(mod.dbp.noOut, d.dbp.noOut)
Random effects are normally distributed. Residuals are homoskedastic. One data point was found to display a cook’s distance > 1 and was removed. Slight deviations from normality appear pre:exp and pre:cont factors. No other outlier is detected. The model is estimated.
an_dbp <- anova(mod.dbp.noOut, ddf = "Kenward-Roger") |> tidy()
an_dbp |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value |
|---|---|---|---|---|---|---|
| timing | 64.932 | 64.932 | 1 | 110.107 | 3.297 | 0.072 |
| condition | 150.254 | 150.254 | 1 | 110.107 | 7.630 | 0.007 |
| timing:condition | 158.279 | 158.279 | 1 | 110.107 | 8.037 | 0.005 |
Models estimated are collated in one table and effects size are computed (\(\omega^2 adjusted\)).
3.1.5 Table models
df_comp <- rbind(an_hr, an_br, an_sbp, an_dbp) |>
mutate(om2 = as_tibble(F_to_omega2(statistic, NumDF, DenDF))[1] |> pull()) |>
mutate(interpretation = as_tibble(interpret_omega_squared(om2))[1] |> pull() )
df_comp$var <- rep(c("HR (BPM)", "BR (CPM)", "SBP (mmHg)", "DBP (mmHg)"),each = 3)
df_comp$var2 <- rep(c("HR", "BR", "SBP", "DBP"),each = 3)
#df_comp$p.corr <- p.adjust(df_comp$p.value, method = "BH")
t_comp <- df_comp |>
select(var, term, sumsq, meansq, NumDF, DenDF, statistic, p.value, om2, interpretation) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
flextable() |>
colformat_double(digits = 3) |>
set_header_labels(var = "Variables", term = "Term", sumsq = "SS",meansq = "MS", NumDF = "Df num", DenDF = "Df den", statistic = "F value", p.value = "p", om2 = "\u03C9\u00B2\u209A", interpretation = "Interpretation") |>
add_header_row(values = c( "Variables", "Term","SS", "MS", "Df num", "Df den", "F value", "p value", "Effect size" ), colwidths = c(1,1,1,1,1,1,1,1,2)
) |>
hline(i = c(3, 6, 9, 12), border = dash_border) |>
merge_v(j = c(1,2,3), part = "body") |>
merge_v(part = "header") |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
autofit() |>
bg(i = ~ p.value < 0.05, j = 8, bg = "yellow") |>
add_footer_lines(values = c("Data from 38 participants. \n Abreviations: HR, Heart rate; BR, Breating rate; SBP, Systolic blood pressure; DBP, Diastolic blood pressure; BPM, Beats per minute; CPM, Cycles per minute; mmHg, millimeter of mercury; SS, Sum of squares; MS, Mean squares; Df, Degrees of freedom."))
# I significant: HR, BR, SBP, DBP, quad, forearm
# main effect significant: duration timing and condition
t_compVariables | Term | SS | MS | Df num | Df den | F value | p value | Effect size | |
|---|---|---|---|---|---|---|---|---|---|
p | ω²ₚ | Interpretation | |||||||
HR (BPM) | timing | 93,307.605 | 93,307.605 | 1 | 111.000 | 2,179.184 | < .001 | 0.951 | large |
condition | 68,467.605 | 68,467.605 | 1 | 111.000 | 1,599.050 | < .001 | 0.934 | large | |
timing:condition | 68,722.526 | 68,722.526 | 1 | 111.000 | 1,605.003 | < .001 | 0.934 | large | |
BR (CPM) | timing | 5,887.605 | 5,887.605 | 1 | 111.000 | 138.522 | < .001 | 0.549 | large |
condition | 3,316.447 | 3,316.447 | 1 | 111.000 | 78.028 | < .001 | 0.405 | large | |
timing:condition | 3,467.605 | 3,467.605 | 1 | 111.000 | 81.585 | < .001 | 0.416 | large | |
SBP (mmHg) | timing | 3,720.421 | 3,720.421 | 1 | 111.000 | 98.744 | < .001 | 0.464 | large |
condition | 3,335.158 | 3,335.158 | 1 | 111.000 | 88.519 | < .001 | 0.436 | large | |
timing:condition | 2,480.237 | 2,480.237 | 1 | 111.000 | 65.828 | < .001 | 0.365 | large | |
DBP (mmHg) | timing | 64.932 | 64.932 | 1 | 110.107 | 3.297 | 0.0721 | 0.020 | small |
condition | 150.254 | 150.254 | 1 | 110.107 | 7.630 | 0.00673 | 0.056 | small | |
timing:condition | 158.279 | 158.279 | 1 | 110.107 | 8.037 | 0.00545 | 0.059 | small | |
Data from 38 participants. | |||||||||
save_as_docx(t_comp, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_comp_modcov.docx")3.1.6 Table pairwise comparisons
fun_emm <- function(mod, var.name){
emm <- emmeans(mod, ~timing*condition)
pairs.t <- pairs(emm, simple = "timing", adjust = "bonferroni") |> as_tibble()
pairs.t.ci <- pairs(emm, simple = "timing", adjust = "bonferroni") |>
confint() |> as_tibble()
pairs.t.full <- full_join(pairs.t, pairs.t.ci)
pairs.c <- pairs(emm, simple = "condition", adjust = "bonferroni") |> as_tibble()
pairs.c.ci <- pairs(emm, simple = "condition", adjust = "bonferroni") |>
confint() |> as_tibble()
pairs.c.full <- full_join(pairs.c, pairs.c.ci)
pairs <- full_join(pairs.t.full, pairs.c.full) |>
select(condition, timing, contrast, estimate, lower.CL, upper.CL, everything())
pairs$var <- rep(var.name, nrow(pairs))
return(pairs)
} # function for pairwise comparisons of significant interactions
l_emm <- list()
l_mod <- list(
mod.hr = mod.hr,
mod.br = mod.br,
mod.sbp = mod.sbp,
mod.dbp = mod.dbp
)
vars <- c("HR (BPM)", "BR (CPM)", "SBP (mmHg)", "DBP (mmHg)")
for (i in 1:length(vars)){
l_emm[[i]] <- fun_emm(l_mod[[i]], vars[i])
}
df_cont <- rbind(l_emm[[1]],l_emm[[2]], l_emm[[3]], l_emm[[4]] )|>
rename(statistic = t.ratio) |>
rowwise() |>
mutate(d_r = as_tibble(t_to_d(statistic, df, paired = T)[1]) |> pull(), interpretation = as.character(as_tibble(interpret_cohens_d(t_to_d(statistic, df, paired = T)))[5] |> pull()))
df_cont$method <- rep("Paired t-test", nrow(df_cont)) # pairwise comparisons of model above are collated in a data frame Normality assessment of of room temperature, room humidity, RPE, SSS, STAI-Y1.
# contrast cont vs exp
d.cont.simple <- data |>
select(ID, condition, RPE, room_temp, room_humidity, SSS, STAI_1) |>
group_by(ID, condition) |>
slice(1) |>
mutate(across(c(RPE, room_temp, room_humidity, SSS, STAI_1), as.numeric)) |>
pivot_longer(cols = c(RPE, room_temp, room_humidity, SSS, STAI_1), names_to = "var", values_to = "values")
# normality assessment
d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
mutate(delta = exp - cont) |>
ggplot(aes(x = delta)) +
geom_density() +
facet_wrap(~var, scales = "free") + my_theme()
d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
mutate(delta = exp - cont) |>
ggplot(aes(sample = delta)) +
geom_qq() +
geom_qq_line() +
facet_wrap(~var, scales = "free") + my_theme()
Room humidity and room temperature values appear normally distributed. Scores of RPE, SSS, and STAI-Y1 violate normaly assumptions.
# RPE
d.rpe <- d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
filter(var == "RPE")
n_pairs <- d.rpe |> # number of non-null pairs
mutate(delta = cont - exp) |>
filter(delta != 0) |>
nrow()
pairs.rpe <- wilcox.test(d.rpe$cont, d.rpe$exp, paired = T, conf.int = T, exact = F) |> tidy() |>
select(estimate, conf.low,conf.high, statistic, p.value) |>
rename(lower.CL = conf.low ) |>
rename(upper.CL = conf.high) |>
mutate()
pairs.rpe$method <- "Wilcoxon signed-rank test"
pairs.rpe$timing <- "post"
pairs.rpe$contrast = "cont - exp"
pairs.rpe$var <- "RPE (6 - 20)"
pairs.rpe$df <- n_pairs
pairs.rpe$d_r <- as_tibble(rank_biserial(d.rpe$cont, d.rpe$exp, paired = T)$r_rank_biserial) |> pull()
pairs.rpe$interpretation <- as.character(as_tibble(interpret_rank_biserial(rank_biserial(d.rpe$cont, d.rpe$exp, paired = T), rules = "gignac2016")$Interpretation) |> pull())
df_cont <- full_join(df_cont, pairs.rpe)
# room_temp
d.rt <- d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
filter(var == "room_temp")
pairs.rt <- t.test(d.rt$cont, d.rt$exp, paired = T) |>
tidy() |>
select(estimate, conf.low, conf.high, statistic, p.value, parameter, method ) |>
rename(df = parameter) |>
rename(lower.CL = conf.low ) |>
rename(upper.CL = conf.high) |>
rowwise() |>
mutate(d_r = as_tibble(t_to_d(statistic, df, paired = T)[1]) |> pull(),
interpretation = as.character(as_tibble(interpret_cohens_d(t_to_d(statistic, df, paired = T)))[5] |> pull()))
pairs.rt$timing <- "pre"
pairs.rt$contrast = "cont - exp"
pairs.rt$var <- "Room temperature (°C)"
df_cont <- full_join(df_cont, pairs.rt)
# room_humidity
d.rh <- d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
filter(var == "room_humidity")
pairs.rh <- t.test(d.rh$cont, d.rh$exp, paired = T) |>
tidy() |> select(estimate, conf.low, conf.high, statistic, p.value, parameter, method ) |>
rename(df = parameter) |>
rename(lower.CL = conf.low ) |>
rename(upper.CL = conf.high)|>
rowwise() |>
mutate(d_r = as_tibble(t_to_d(statistic, df, paired = T)[1]) |> pull(),
interpretation = as.character(as_tibble(interpret_cohens_d(t_to_d(statistic, df, paired = T)))[5] |> pull()))
pairs.rh$timing <- "pre"
pairs.rh$contrast = "cont - exp"
pairs.rh$var <- "Room humidity (%)"
df_cont <- full_join(df_cont, pairs.rh)
# STAI1
d.stai1 <- d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
filter(var == "STAI_1")
n_pairs <- d.stai1 |> # number of non-null pairs
mutate(delta = cont - exp) |>
filter(delta != 0) |>
nrow()
pairs.stai1 <- wilcox.test(d.stai1$cont, d.stai1$exp, paired = T, conf.int = T, exact = F) |> tidy() |>
select(estimate, conf.low,conf.high, statistic, p.value) |>
rename(lower.CL = conf.low ) |>
rename(upper.CL = conf.high) |>
mutate()
pairs.stai1$method <- "Wilcoxon signed-rank test"
pairs.stai1$timing <- "pre"
pairs.stai1$contrast = "cont - exp"
pairs.stai1$var <- "STAI-Y1 (/80)"
pairs.stai1$df <- n_pairs
pairs.stai1$d_r <- as_tibble(rank_biserial(d.stai1$cont, d.stai1$exp, paired = T)$r_rank_biserial) |> pull()
pairs.stai1$interpretation <- as.character(as_tibble(interpret_rank_biserial(rank_biserial(d.stai1$cont, d.stai1$exp, paired = T), rules = "gignac2016")$Interpretation) |> pull())
df_cont <- full_join(df_cont, pairs.stai1)
# SSS
d.sss <- d.cont.simple |>
pivot_wider(names_from = condition, values_from = values) |>
filter(var == "SSS")
n_pairs <- d.sss |> # number of non-null pairs
mutate(delta = cont - exp) |>
filter(delta != 0) |>
nrow()
pairs.sss <- wilcox.test(d.sss$cont, d.sss$exp, paired = T, conf.int = T, exact = F) |> tidy() |>
select(estimate, conf.low,conf.high, statistic, p.value) |>
rename(lower.CL = conf.low ) |>
rename(upper.CL = conf.high) |>
mutate()
pairs.sss$method <- "Wilcoxon signed-rank test"
pairs.sss$timing <- "pre"
pairs.sss$contrast = "cont - exp"
pairs.sss$var <- "SSS (/7)"
pairs.sss$df <- n_pairs
pairs.sss$d_r <- as_tibble(rank_biserial(d.sss$cont, d.sss$exp, paired = T)$r_rank_biserial) |> pull()
pairs.sss$interpretation <- as.character(as_tibble(interpret_rank_biserial(rank_biserial(d.sss$cont, d.sss$exp, paired = T), rules = "gignac2016")$Interpretation) |> pull())
df_cont <- full_join(df_cont, pairs.sss) Pairwise comparisons of significant interactions or main effects from models and paired tests are collated in one table.
t1_contrasts <- df_cont |>
select(var, timing, condition, contrast, estimate,lower.CL, upper.CL, SE,method,df, statistic,p.value, d_r, interpretation ) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
flextable() |>
set_header_labels(var = "Variables", timing = "Timing", condition = "Condition", contrast = "Contrast", estimate = "Difference", lower.CL = "low", upper.CL = "high", SE = "SE", method = "Method", df = "Df", statistic = "Statistic", p.value = "p", d_r = "d~z~ / r~rb~", interpretation = "Interpretation") |>
colformat_md(part = "header", j = 13) |>
add_header_row(values = c("Variables", "Timing", "Condition", "Contrast", "Difference", "CI 95%", "SE", "Method", "Df", "Statistic", "p value", "Effect sizes" ), colwidths = c(1,1,1,1,1,2,1,1,1,1,1, 2)) |>
colformat_double(digits = 3) |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
merge_v(j = c(1), part = "body") |>
bg(i = ~ p.value < 0.05, j = 12, bg = "yellow") |>
merge_v(part = "header") |>
hline(i = c(4, 8, 12, 16, 17,18,19,20), border = dash_border) |>
autofit()|>
add_footer_lines(values = c("Data from 38 participants. \n Method refers to the statistical test used for the contrast. Degrees of freedeom in case of wilcoxon signed rank test refers to number of non-null pairs. \n Abreviations: HR, Heart rate; BR, Breating rate; SBP, Systolic blood pressure; DBP, Diastolic blood pressure; T, temperature; RPE, Rate of perceived exertion; STAI-Y1, State anxiety questionnaire; SSS, Stanford sleepiness scale; BPM, Beats per minute; CPM, Cycles per minute; mmHg, millimeter of mercury; °C, Celcius degrees; SE, Standard error; CI 95%, 95% confidence interval; Df, Degrees of freedom; dz, Cohen'd for paired samples; rrb, rank-biserial correlation."))
t1_contrasts Variables | Timing | Condition | Contrast | Difference | CI 95% | SE | Method | Df | Statistic | p value | Effect sizes | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
low | high | p | dz / rrb | Interpretation | |||||||||
HR (BPM) | cont | pre - post | -7.026 | -10.001 | -4.052 | 1.501 | Paired t-test | 111.000 | -4.681 | < .001 | -0.444 | small | |
exp | pre - post | -92.079 | -95.054 | -89.104 | 1.501 | Paired t-test | 111.000 | -61.337 | < .001 | -5.822 | large | ||
pre | cont - exp | 0.079 | -2.896 | 3.054 | 1.501 | Paired t-test | 111.000 | 0.053 | 0.958 | 0.005 | very small | ||
post | cont - exp | -84.974 | -87.948 | -81.999 | 1.501 | Paired t-test | 111.000 | -56.604 | < .001 | -5.373 | large | ||
BR (CPM) | cont | pre - post | -2.895 | -5.858 | 0.069 | 1.496 | Paired t-test | 111.000 | -1.935 | 0.0555 | -0.184 | very small | |
exp | pre - post | -22.000 | -24.964 | -19.036 | 1.496 | Paired t-test | 111.000 | -14.709 | < .001 | -1.396 | large | ||
pre | cont - exp | 0.211 | -2.753 | 3.174 | 1.496 | Paired t-test | 111.000 | 0.141 | 0.888 | 0.013 | very small | ||
post | cont - exp | -18.895 | -21.858 | -15.931 | 1.496 | Paired t-test | 111.000 | -12.633 | < .001 | -1.199 | large | ||
SBP (mmHg) | cont | pre - post | -1.816 | -4.606 | 0.975 | 1.408 | Paired t-test | 111.000 | -1.289 | 0.2 | -0.122 | very small | |
exp | pre - post | -17.974 | -20.764 | -15.183 | 1.408 | Paired t-test | 111.000 | -12.764 | < .001 | -1.211 | large | ||
pre | cont - exp | -1.289 | -4.080 | 1.501 | 1.408 | Paired t-test | 111.000 | -0.916 | 0.362 | -0.087 | very small | ||
post | cont - exp | -17.447 | -20.238 | -14.657 | 1.408 | Paired t-test | 111.000 | -12.390 | < .001 | -1.176 | large | ||
DBP (mmHg) | cont | pre - post | 0.737 | -1.440 | 2.914 | 1.099 | Paired t-test | 111.000 | 0.671 | 0.504 | 0.064 | very small | |
exp | pre - post | -3.947 | -6.125 | -1.770 | 1.099 | Paired t-test | 111.000 | -3.593 | < .001 | -0.341 | small | ||
pre | cont - exp | 0.053 | -2.125 | 2.230 | 1.099 | Paired t-test | 111.000 | 0.048 | 0.962 | 0.005 | very small | ||
post | cont - exp | -4.632 | -6.809 | -2.454 | 1.099 | Paired t-test | 111.000 | -4.215 | < .001 | -0.400 | small | ||
RPE (6 - 20) | post | cont - exp | -8.500 | -9.000 | -8.000 | Wilcoxon signed-rank test | 38.000 | 0.000 | < .001 | -1.000 | large | ||
Room temperature (°C) | pre | cont - exp | 0.516 | -0.042 | 1.074 | Paired t-test | 37.000 | 1.873 | 0.069 | 0.308 | small | ||
Room humidity (%) | pre | cont - exp | 1.632 | -0.556 | 3.819 | Paired t-test | 37.000 | 1.511 | 0.139 | 0.248 | small | ||
STAI-Y1 (/80) | pre | cont - exp | -2.000 | -5.000 | 1.000 | Wilcoxon signed-rank test | 31.000 | 176.000 | 0.161 | -0.290 | moderate | ||
SSS (/7) | pre | cont - exp | -0.000 | -0.500 | 1.000 | Wilcoxon signed-rank test | 24.000 | 148.500 | 0.975 | -0.010 | very small | ||
Data from 38 participants. | |||||||||||||
save_as_docx(t1_contrasts, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_mcov_contrasts.docx")
#print(t1_contrasts, preview = "docx")3.1.7 Plots
Data used in model is plotted along with significance results.
# HR
d.hr.sum <- d.mod.cov(hr)$dsum
plot.hr <- plot_prepost_an(d.hr, t = "Heart rate", u = "BPM", df_an = df_comp, v = "HR", y = "value") +
showSignificance( c(1.05,2.1), max(d.hr.sum$cih) + 0.02*max(d.hr.sum$cih), "*", width = -0.01, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( c(0.9,1.95), tn(d.hr.sum$cih, 2) + 0.05*tn(d.hr.sum$cih,2), "*", width = -0.01, segmentParams = list(col = "#00468BFF"), textParams = list(col = "#00468BFF")) +
showSignificance( 2.2, c(min(d.hr.sum$cil),max(d.hr.sum$cih)), "*", width = -0.01) +
scale_y_continuous(limits = c(20, 200), breaks = c(40, 80, 120, 160)) +
scale_x_discrete(breaks = "") +
labs(x = "")
# BR
d.br.sum <- d.mod.cov(br)$dsum
plot.br <- plot_prepost_an(d.br, t = "Breating rate", u = "CPM", df_an = df_comp, v = "BR", y = "value")+ showSignificance( c(1.05,2.1), max(d.br.sum$cih) + 0.05*max(d.br.sum$cih), "*", width = -0.01, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( 2.2, c(min(d.br.sum$cil),max(d.br.sum$cih)), "*", width = -0.01) +
scale_x_discrete(breaks = "") +
labs(x = "")
# SBP
d.sbp.sum <- d.mod.cov(sbp)$dsum
plot.sbp <- plot_prepost_an(d.sbp, t = "Systolic blood pressure", u = "mmHg", df_an = df_comp, v = "SBP", y = "value") + showSignificance( c(1.05,2.1), max(d.sbp.sum$cih) + 0.05*max(d.sbp.sum$cih), "*", width = -0.01, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( 2.2, c(min(d.sbp.sum$cil),max(d.sbp.sum$cih)), "*", width = -0.01)
#DBP
d.dbp.sum <- d.mod.cov(dbp)$dsum
plot.dbp <- plot_prepost_an(d.dbp, t = "Diastolic blood pressure", u = "mmHg", df_an = df_comp, v = "DBP", y = "value")+ showSignificance( c(1.05,2.1), max(d.dbp.sum$cih) + 0.05*max(d.dbp.sum$cih), "*", width = -0.01, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( 2.2, c(min(d.dbp.sum$cil),max(d.dbp.sum$cih)), "*", width = -0.01)
p.cov <- ggarrange(plot.hr, plot.br, plot.sbp, plot.dbp, ncol = 2, nrow = 2, common.legend = T, legend = "top", align = "v", labels = "AUTO" )
ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/pcov.png", plot = p.cov, height = 7, width = 7,dpi = 600)
p.cov
3.1.8 Table results summary
Described results are summarized in a table along with summary statistics.
d.covar.sum <- d_cov_long |>
select(timing, condition, hr, br, sbp, dbp, ID) |>
pivot_longer(cols = c(hr, br, sbp, dbp),
names_to = c("var"),values_to = "values" ) |>
group_by(var, condition, timing) |>
summarise(m = mean(values),
sd = sd(values),
md = median(values),
iqr = IQR(values))
d.cont.simple.sum <- d.cont.simple |>
group_by(var, condition) |>
summarise(m = mean(values),
sd = sd(values),
md = median(values),
iqr = IQR(values))
df_t2 <- full_join(d.covar.sum,d.cont.simple.sum ) |>
mutate(across(c(m, sd), round, 2)) |>
mutate(
msd = glue("{m} ± {sd}"),
med_iqr = glue("{md} [{iqr}]^a^")) |>
select(-m, -sd, -md, -iqr) |>
pivot_wider(names_from = condition, values_from = c(msd, med_iqr)) |>
mutate(msd_exp = if_else(var %in% c("RPE", "SSS", "STAI_1"), med_iqr_exp, msd_exp),
msd_cont = if_else(var %in% c("RPE", "SSS", "STAI_1"), med_iqr_cont, msd_cont)
) |>
rename(cont = msd_cont) |>
rename(exp = msd_exp) |>
select(-starts_with("med")) |>
mutate(timing = if_else(var %in% c("room_humidity", "room_temp", "SSS", "STAI_1"), "pre", timing),
timing = if_else(var == "RPE", "post", timing)) |>
pivot_wider(names_from = timing, values_from = c(exp, cont)) |>
select(var, cont_pre, exp_pre, cont_post, exp_post) |>
mutate(var = factor(var, levels = c("hr", "br", "sbp", "dbp", "RPE", "SSS", "STAI_1", "room_humidity", "room_temp"), labels = c("HR (BPM)", "BR (CPM)", "SBP (mmHg)", "DBP (mmHg)", "RPE (6 - 20)", "SSS (/7)", "STAI-Y1 (/80)", "Room humidity (%)","Room temperature (°C)"))) |>
arrange(var)
p_inter <- df_comp |> filter(term == "timing:condition") |>
mutate(across(c(NumDF, DenDF, statistic), round, 2)) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
mutate(res = glue("F~{NumDF},{DenDF}~ = {statistic}")) |>
mutate(var = factor(var, levels= c("HR (BPM)", "BR (CPM)", "SBP (mmHg)", "DBP (mmHg)"))) |>
arrange(var)
p_cont <- df_cont |>
filter(var %in% c("RPE (6 - 20)", "SSS (/7)", "STAI-Y1 (/80)", "Room temperature (°C)", "Room humidity (%)")) |>
mutate(var = factor(var, levels = c("RPE (6 - 20)", "SSS (/7)", "STAI-Y1 (/80)", "Room temperature (°C)", "Room humidity (%)"))) |>
mutate(across(c(statistic, df), round, 2)) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
mutate(stat = if_else(method == "Paired t-test", "t", "V")) |>
mutate(res = glue("{stat}~{df}~ = {statistic}")) |>
arrange(var)
df_t2$stat <- c(p_inter$res, p_cont$res)
df_t2$p.value <- c(p_inter$p.value, p_cont$p.value)
t2 <- df_t2 |>
flextable() |>
set_header_labels(var = "Variables", cont_pre = "Control", exp_pre = "Exercise", cont_post = "Control", exp_post = "Exercise", stat = "Stat~df~", p.value = "p") |>
add_header_row(values = c("Variables", "Pre", "Post", "Timing x Condition interaction / Pairwise comparison"), colwidths = c(1, 2,2,2)) |>
bold(part = "header") |>
merge_v(part = "header") |>
align(part = "header", align = "center") |>
vline(j = c(1,3,5), border = dash_border) |>
autofit() |>
colformat_md(part = "all", j = 1:7) |>
add_footer_lines(values = c("a Median [interquartile range]. Data from 38 participants. \n Abreviations: HR, Heart rate; BR, Breating rate; SBP, Systolic blood pressure; DBP, Diastolic blood pressure; T, temperature; BPM, Beats per minute; CPM, Cycles per minute; mmHg, millimeter of mercury; RPE, Rate of perceived exertion; °C, Celsius degrees; Df, Degrees of freedom"))
t2Variables | Pre | Post | Timing x Condition interaction / Pairwise comparison | |||
|---|---|---|---|---|---|---|
Control | Exercise | Control | Exercise | Statdf | p | |
HR (BPM) | 61.58 ± 8.87 | 61.5 ± 11.12 | 68.61 ± 9.71 | 153.58 ± 9.2 | F1,111 = 1605 | < .001 |
BR (CPM) | 13.32 ± 3.63 | 13.11 ± 3.09 | 16.21 ± 4.79 | 35.11 ± 12.04 | F1,111 = 81.58 | < .001 |
SBP (mmHg) | 113.45 ± 9.84 | 114.74 ± 7.1 | 115.26 ± 8.08 | 132.71 ± 10.52 | F1,111 = 65.83 | < .001 |
DBP (mmHg) | 65.32 ± 5.71 | 65.26 ± 5.86 | 64.58 ± 6.82 | 69.21 ± 8.19 | F1,110.11 = 8.04 | 0.00545 |
RPE (6 - 20) | 6 [1]a | 15 [2]a | V38 = 0 | < .001 | ||
SSS (/7) | 2 [1]a | 2 [1]a | V24 = 148.5 | 0.975 | ||
STAI-Y1 (/80) | 29.5 [9]a | 29 [10.75]a | V31 = 176 | 0.161 | ||
Room humidity (%) | 47.16 ± 7.79 | 45.53 ± 7.52 | t37 = 1.87 | 0.069 | ||
Room temperature (°C) | 21.8 ± 1.52 | 21.28 ± 1.61 | t37 = 1.51 | 0.139 | ||
a Median [interquartile range]. Data from 38 participants. | ||||||
#print(t1, preview = "docx")
save_as_docx(t2, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t1.docx")3.2 Carry over effect
Before to compute the model, we test for potential carry-over effects. Although unlikely, it would be possible that the carry over effect of receiving exercise in the first session impacts more the results on the second session than the other way around.
First, we visualise the absolute change (post intervention - pre intervention) across sequences (exercise/control vs control/exercise) and conditions (exercise vs control)
# inclusion of sequence factor in data set
# computation of mean absolute change
s_ppt_q <- data_long |>
mutate(seq = glue("{session}_{condition}")) |>
rowwise() |>
mutate(test_value = if_else(test_mod %in% c("ppt", "cdt"), logt(test_value), test_value)) |>
select(ID, order, test_mod,timing, site, session , condition, test_value, seq) |>
mutate(seq = factor(seq)) |>
mutate(seq_cont_exp = if_else(seq %in% c("session_1_cont", "session_2_exp"), "cont_exp", "exp_cont")) |>
select(-seq) |>
mutate(seq = seq_cont_exp) |>
select(!seq_cont_exp)|>
group_by(ID, condition, test_mod, site, timing) |>
mutate(tm = mean(test_value)) |>
slice(1) |>
select(-test_value) |>
ungroup() |>
pivot_wider(
names_from = timing,
values_from = tm
) |>
mutate(abc = post - pre) |>
ungroup() |>
group_by(session, seq, condition, test_mod, site) |>
summarize(
m = mean(abc, na.rm = T),
s = sd(abc, na.rm = T)
)
co_seq <- s_ppt_q |> # sequence plot
ggplot(aes(x = session, y = m, col = seq, group = seq)) +
geom_point() +
geom_line() +
facet_grid(rows = vars(test_mod), cols = vars(site), scale = "free_y") +
labs(title = "Absolute change ~ session * sequence")
co_cond <- s_ppt_q |> # condition plot
ggplot(aes(x = session, y = m, col = condition, group = condition)) +
geom_point() +
geom_line() +
facet_grid(rows = vars(test_mod), cols = vars(site), scale = "free_y") +
labs(title = "Absolute change ~ session * condition")
co_seq
co_cond
Mean absolute change values appear greater in exp_cont vs cont_exp in both sessions for hpt quadriceps, and pp quadriceps and forearm; they appear greater in cont_exp vs exp_cont for ppt forearm. This suggests a possible carry-over effects of sequences. A slight increase of mean absolute change values in observed across sessions for ppt quad, ppt forearm, cdt quad, cdt forearm.
Possible carry over effects are tested by comparing the sum of absolute change in each sequence using t test for each test modality:site. P values are corrected for false discovery rate.
# this is a script to understand how to take into account carry over effects in within subj design
# it is based on Shen et al., Estimate Carryover Effect in Clinical Trial Crossover Designs
# Vladimir Aron
# 20240521
carry_test <- function(data,s , test_m ) {
## data <- data_long
# s <- "no site"
# test_m <- "adt"
co <- data |>
rowwise() |>
mutate(test_value = if_else(test_mod %in% c("ppt", "cdt"), logt(test_value), test_value)) |>
mutate(seq = glue("{session}_{condition}")) |>
mutate(seq = factor(seq)) |>
mutate(seq_cont_exp = if_else(seq %in% c("session_1_cont", "session_2_exp"), "cont_exp", "exp_cont")) |>
select(-seq) |>
mutate(seq = seq_cont_exp) |>
select(!seq_cont_exp)
# filter(site == "quad", test_mod == "ppt")
qppt <- co |>
select(ID, order, test_mod,timing, site, session , condition, test_value, seq) |>
group_by(ID, condition, test_mod, site, timing) |>
mutate(tm = mean(test_value)) |>
slice(1) |>
select(!test_value) |>
ungroup() |>
filter(site == s, test_mod == test_m) |>
ungroup() |>
pivot_wider(
names_from = timing,
values_from = tm
) |>
mutate(abc = post - pre) |>
ungroup() |>
select(ID, seq, session, condition, abc) |>
arrange(ID)
qppt.s <- qppt |>
group_by(ID) |>
mutate(sum = sum(abc))
qppt.s.ab <- subset(qppt.s, seq == "exp_cont") |>
group_by(ID) |>
slice(1)
qppt.s.ba <- subset(qppt.s, seq == "cont_exp") |>
group_by(ID) |>
slice(1)
# returns a tibble
t_res <- t.test(qppt.s.ab$sum, qppt.s.ba$sum)
res <- broom::tidy(t_res)
}df <- data.frame(
site = c(rep("quad", 4), rep("forearm", 4), "no site"),
test_mod = c(rep(levels(data_long$test_mod)[-5],2),"adt")
)
df_e <- rbind(
carry_test(data_long, s = df[1,1], test_m = df[1,2]),
carry_test(data_long,s = df[1,1], test_m= df[2,2]),
carry_test(data_long,s = df[1,1], test_m= df[3,2]),
carry_test(data_long,s = df[1,1], test_m = df[4,2]),
carry_test(data_long,s = df[5,1], test_m = df[1,2]),
carry_test(data_long,s = df[5,1], test_m = df[2,2]),
carry_test(data_long,s = df[5,1], test_m = df[3,2]),
carry_test(data_long,s = df[5,1], test_m = df[4,2]),
carry_test(data_long,s = df[9,1], test_m = df[9,2])
)
co_df <- cbind(df, df_e) |>
rename("exp_cont" = "estimate1",
"cont_exp" = "estimate2") |>
select(-parameter, -method, -alternative, conf.low, conf.high, -statistic) |> ungroup()
co_df$p_corr <- p.adjust(co_df$p.value, method = "BH")
t_carry <- co_df |>
mutate(test_mod = toupper(test_mod)) |>
select(site, test_mod, exp_cont, cont_exp, estimate, conf.low, conf.high, p.value, p_corr) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value),
p_corr = pvalue_format_table(p_corr)) |>
mutate(site = factor(site, levels = c("forearm", "quad", "nosite") , labels = c("Forearm", "Quadriceps", "No site"))) |>
flextable() |>
set_header_labels(site = "Site", test_mod = "Test modalities", exp_cont = "Exp/Cont", cont_exp = "Cont/Exp", estimate = "Sum difference", conf.low = "low", conf.high = "high", p.value = "p", p_corr = "p corr") |>
add_header_row(
values = c("Site", "Test modalities","Sum absolute change", "Sum difference","CI 95%", "p", "p corr"),
colwidths = c(1, 1, 2,1,2, 1, 1)
) %>%
colformat_double(digits = 3) |>
bold(part = "header") |>
merge_at(j = 1, i = 1:4) %>%
merge_at(j = 1, i = 5:8,part = "body") %>%
merge_v(part = "header") |>
hline(i = c(4,8), border = dash_border) %>%
autofit() |>
add_footer_lines(values = c("Data from 38 participants.\n Units: CDT, log; HPT, °C; PP, VAS 100 - 0; PPT, log.\n Abbreviations: CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Mechanical pinprick rating; PPT, Pressure pain threshold; CI 95%, 95% confidence interval; p corr, p value ajdusted with false discovery rate."))
ts2 <- set_table_properties(t_carry, layout = "autofit")
#print(ts2, preview = "docx")
ts2Site | Test modalities | Sum absolute change | Sum difference | CI 95% | p | p corr | ||
|---|---|---|---|---|---|---|---|---|
Exp/Cont | Cont/Exp | low | high | |||||
Quadriceps | PPT | 0.001 | 0.007 | -0.006 | -0.069 | 0.057 | 0.848 | 0.968 |
PP | -1.026 | -5.330 | 4.304 | -4.702 | 13.311 | 0.339 | 0.968 | |
HPT | -0.056 | -0.758 | 0.702 | -0.821 | 2.225 | 0.356 | 0.968 | |
CDT | -0.036 | -0.031 | -0.005 | -0.129 | 0.119 | 0.938 | 0.968 | |
Forearm | PPT | -0.013 | 0.011 | -0.024 | -0.118 | 0.070 | 0.607 | 0.968 |
PP | 2.147 | -9.073 | 11.220 | 0.507 | 21.932 | 0.0406 | 0.365 | |
HPT | 0.282 | -0.158 | 0.440 | -0.962 | 1.842 | 0.528 | 0.968 | |
CDT | 0.027 | 0.025 | 0.002 | -0.103 | 0.107 | 0.968 | 0.968 | |
ADT | 0.001 | -0.001 | 0.002 | -0.017 | 0.022 | 0.821 | 0.968 | |
Data from 38 participants. | ||||||||
save_as_docx(ts2, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_carry.docx")No significant carry-over effect is observed.
3.3 Reliability of pre values
Reliability of pre intervention values for each test modality and site is assessed by computing ICC(3,k)
# function to compute reliability of pre values
# 20240527
# Vladimir Aron
rel_pre <- function(data, anov = T) {
# function that analyze reliaiblity of wp1b EIH data using ICC function from psych package
# outputs tables with ICC
# input: data set (data_long) and anov relates to computation method in ICC function: anova vs lme; default: ANOVA
data.rel <- function(data, s, tm) {
# mean values of three measurements are computed
d <- data |>
filter(timing == "pre") |>
group_by(ID, session, test_mod, site) |>
mutate(m = mean(test_value, na.rm = T)) |>
rowwise() |>
mutate(m = if_else(test_mod %in% c("cdt", "ppt"), logt(m), m))
data.mod.s1 <- d |>
filter(site == s, test_mod == tm, session == "session_1") |>
group_by(ID) |>
slice(1)
data.mod.s2 <- d |>
filter(site == s, test_mod == tm, session == "session_2") |>
group_by(ID) |>
slice(1)
d <- cbind(data.mod.s1$m, data.mod.s2$m) # returns a df of individual values for each session
return(d)
}
data.rel.qppt <- data.rel(data, s = "quad", tm ="ppt")
data.rel.qcdt <- data.rel(data, s = "quad", tm ="cdt")
data.rel.qhpt <- data.rel(data, s = "quad", tm ="hpt")
data.rel.qpp <- data.rel(data, s = "quad", tm ="pp")
data.rel.fppt <- data.rel(data, s = "forearm", tm ="ppt")
data.rel.fcdt <- data.rel(data, s = "forearm", tm ="cdt")
data.rel.fhpt <- data.rel(data, s = "forearm", tm ="hpt")
data.rel.fpp <- data.rel(data, s = "forearm", tm ="pp")
data.rel.adt <- data.rel(data, s = "no site", tm = "adt")
if (anov == T) {
icc.res.qppt <- psych::ICC(data.rel.qppt, lmer = F)
icc.res.qcdt <- psych::ICC(data.rel.qcdt, lmer = F)
icc.res.qhpt <- psych::ICC(data.rel.qhpt, lmer = F)
icc.res.qpp <- psych::ICC(data.rel.qpp, lmer = F)
icc.res.fppt <- psych::ICC(data.rel.fppt, lmer = F)
icc.res.fcdt <- psych::ICC(data.rel.fcdt, lmer = F)
icc.res.fhpt <- psych::ICC(data.rel.fhpt, lmer = F)
icc.res.fpp <- psych::ICC(data.rel.fpp, lmer = F)
icc.res.adt <- psych::ICC(data.rel.adt, lmer = F)
comput <- "anova"
} else {
icc.res.qppt <- psych::ICC(data.rel.qppt, lmer = T)
icc.res.qcdt <- psych::ICC(data.rel.qcdt, lmer = T)
icc.res.qhpt <- psych::ICC(data.rel.qhpt, lmer = T)
icc.res.qpp <- psych::ICC(data.rel.qpp, lmer = T)
icc.res.fppt <- psych::ICC(data.rel.fppt, lmer = T)
icc.res.fcdt <- psych::ICC(data.rel.fcdt, lmer = T)
icc.res.fhpt <- psych::ICC(data.rel.fhpt, lmer = T)
icc.res.fpp <- psych::ICC(data.rel.fpp, lmer = T)
icc.res.adt <- psych::ICC(data.rel.adt, lmer = T)
comput <- "lme"
}
icc_df <- rbind(
icc.res.qppt$results |> filter(type == "ICC3k"),
icc.res.qcdt$results |> filter(type == "ICC3k"),
icc.res.qhpt$results |> filter(type == "ICC3k"),
icc.res.qpp$results |> filter(type == "ICC3k"),
icc.res.fppt$results |> filter(type == "ICC3k"),
icc.res.fcdt$results |> filter(type == "ICC3k"),
icc.res.fhpt$results |> filter(type == "ICC3k"),
icc.res.fpp$results |> filter(type == "ICC3k"),
icc.res.adt$results |> filter(type == "ICC3k")
)
icc_df$site <- c(rep("quad",4), rep("forearm",4), "nosite")
icc_df$test_mod <- c(rep(c("PPT", "CDT", "HPT", "PP"),2), "ADT")
row.names(icc_df) <- NULL
t.icc <- icc_df |>
select(site,test_mod, ICC, `lower bound`, `upper bound`) |>
mutate(across(is.numeric, round, 2)) |>
kbl(caption = paste("ICC values (", comput, ") at all sites & test modalities"),
col.names = c(site = "Site", test_mod = "Test modality", ICC = "ICC(3C,k)", `lower bound` = "CI low", `upper bound` = "CI high")) |>
kable_minimal() |>
collapse_rows(columns = 1, valign = "top")
results <- list(
"t.icc" = t.icc,
"icc.df" = icc_df
)
return(results)
}t_rel <- rel_pre(data_long)$icc.df |>
select(site,test_mod, ICC, `lower bound`, `upper bound`) |>
mutate(site = factor(site, levels = c("forearm", "quad", "nosite"), labels = c("Forearm", "Quadriceps", "No site"))) |>
flextable() |>
colformat_double(digits = 2) |>
set_header_labels(site = "Site", test_mod = "Test modalities", ICC = "ICC" ,`lower bound` = "low",`upper bound` = "high") |>
add_header_row(values = c("Site", "Test modalities", "ICC", "CI 95%"), colwidths = c(1,1,1,2))|>
bold(part = "header") |>
merge_at(j = 1, i = 1:4) %>%
merge_at(j = 1, i = 5:8,part = "body") %>%
merge_v(part = "header") |>
hline(i = c(4,8), border = dash_border) %>%
autofit() |>
add_footer_lines(values = c("Data from 38 participants.\n Abbreviations: CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Mechanical pinprick rating; PPT, Pressure pain threshold; ICC, Intraclass correlation coefficient; CI 95%, 95% confidence interval."))
ts3 <- set_table_properties(t_rel, layout = "autofit")
#print(ts3, preview = "docx")
ts3Site | Test modalities | ICC | CI 95% | |
|---|---|---|---|---|
low | high | |||
Quadriceps | PPT | 0.88 | 0.77 | 0.94 |
CDT | 0.84 | 0.70 | 0.92 | |
HPT | 0.95 | 0.90 | 0.97 | |
PP | 0.86 | 0.73 | 0.93 | |
Forearm | PPT | 0.84 | 0.69 | 0.92 |
CDT | 0.88 | 0.78 | 0.94 | |
HPT | 0.83 | 0.67 | 0.91 | |
PP | 0.91 | 0.82 | 0.95 | |
No site | ADT | 0.84 | 0.70 | 0.92 |
Data from 38 participants. | ||||
save_as_docx(ts3, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_rel.docx")All ICC are > 0.7 suggesting that sensory testing is reliable.
3.4 Skin temperature
Finally, we analyse the changes in skin temperature before and after each condition at both forearms and quadriceps.
# long format skin temperature data set
d_long_temp <- data |>
select(starts_with("temp"), condition, ID) |>
select(-temp_room) |>
mutate(across(starts_with("temp"), as.numeric)) |>
pivot_longer(cols = starts_with("temp"),
names_to = c("temp_t","site_t", "time_t","dominance_t"),
names_pattern = "(^\\w+)_(\\w+)_(\\w+)_(\\w+)",
values_to = c("temp")) |>
mutate(site = paste0(site_t, "_", dominance_t)) |>
select(- site_t, -dominance_t, -temp_t) |>
mutate(
time_t = factor(time_t, levels = c("pre1", "pre2", "post1", "post2")),
site = factor(site),
condition = factor(condition, levels = c("exp", "cont"))) |>
group_by(ID, condition, site, time_t) |>
slice(1) |>
na.omit()
# mean post - pre temperature values
delta_temp <- d_long_temp |>
pivot_wider(names_from = "time_t", values_from = "temp") |>
mutate(pre = (pre1 + pre2)/2,
post = (post1 + post2)/2,
abc_temp = post - pre) |>
mutate(condition = factor(condition, levels = c("cont", "exp"), labels = c("Control", "Exercise")),
site = factor(site, levels = levels(factor(d_long_temp$site)), labels = c("Dominant forearm", "Non-dominant forearm", "Dominant quadriceps", "Non-dominant quadriceps"))
)
# long format mean temperature
d_temp_mean <- delta_temp |>
pivot_longer(cols = c("pre", "post"), names_to = "timing", values_to = "temp") |>
mutate(timing = factor(timing, levels = c("pre", "post"), labels = c("Pre", "Post")))
# delta temperature summary data set
delta_temp_sum <- delta_temp |>
group_by(condition, site) |>
summarise(m = mean(abc_temp),
sd = sd(abc_temp),
cil = ci_low(abc_temp),
cih = ci_high(abc_temp)) |>
mutate(cond_s = paste(condition, site)) |>
arrange(cond_s) %>%
mutate(unique_cond = paste(condition, site)) %>%
arrange(unique_cond)
# mean temp summary data set
d_temp_sum <- d_temp_mean |>
group_by(site, condition, timing) |>
summarise(
m = mean(temp),
sd = sd(temp),
cil = ci_low(temp),
cih = ci_high(temp)
)
# pre value normality
delta_temp |>
select(condition, site, pre) |>
ungroup() |>
group_by(condition, site) |>
shapiro_test(pre)delta_temp |>
ggplot(aes(sample = pre)) +
stat_qq() +
stat_qq_line(col = "red") +
facet_grid(cols = vars(condition), rows = vars(site)) +
my_theme()
delta_temp |>
ggplot(aes(x= pre)) +
geom_histogram(aes(y = after_stat(density))) +
geom_density(fill = "red", alpha = 0.5) +
facet_grid(cols = vars(condition), rows = vars(site)) +
my_theme()
# post value normality
delta_temp |>
select(condition, site, post) |>
ungroup() |>
group_by(condition, site) |>
shapiro_test(post)delta_temp |>
ggplot(aes(sample = post)) +
stat_qq() +
stat_qq_line(col = "red") +
facet_grid(cols = vars(condition), rows = vars(site)) +
my_theme()
delta_temp |>
ggplot(aes(x= post)) +
geom_histogram(aes(y = after_stat(density))) +
geom_density(fill = "red", alpha = 0.5) +
facet_grid(cols = vars(condition), rows = vars(site)) +
my_theme()
# absolute change normality
delta_temp |>
select(condition, site, abc_temp) |>
ungroup() |>
group_by(condition, site) |>
shapiro_test(abc_temp)delta_temp |>
ggplot(aes(sample = abc_temp)) +
stat_qq() +
stat_qq_line(col = "red") +
facet_grid(cols = vars(condition), rows = vars(site)) +
my_theme()
delta_temp |>
ggplot(aes(x= abc_temp)) +
geom_histogram(aes(y = after_stat(density))) +
geom_density(fill = "red", alpha = 0.5) +
facet_grid(cols = vars(condition), rows = vars(site)) +
my_theme()
All levels site:condition appear normally distributed except for dominant_quadriceps:control where are negative skew is observed, potentially due to an outlier value. The outlier value is kept. Summary statistics are computed:
t_summary_temp <- d_long_temp |>
pivot_wider(names_from = "time_t", values_from = "temp") |>
mutate(pre = (pre1 + pre2)/2,
post = (post1 + post2)/2,
abc_temp = post - pre) |>
select(-pre1, -pre2, -post1, -post2) |>
pivot_longer(cols = c("pre", "post", "abc_temp"), names_to = "timing", values_to = "temp") |>
mutate(cond_timing = paste0(timing,"_", condition)) |>
select(-condition, -timing) |>
group_by(site, cond_timing) |>
summarise(m = specify_decimal(mean(temp, na.rm = T), 2),
sd = specify_decimal(sd(temp, na.rm = T), 2),
med = specify_decimal(median(temp, na.rm = T), 2),
iqr = specify_decimal(IQR(temp, na.rm = T), 2)
) |>
mutate(msd = glue("{m} ± {sd}")
) |>
select(-m, -sd, - med, -iqr) |>
pivot_wider(names_from = cond_timing, values_from = msd) |>
select(site, pre_cont, post_cont, abc_temp_cont, pre_exp, post_exp, abc_temp_exp) |>
mutate(
site = factor(site, levels = levels(factor(d_long_temp$site)), labels = c("Dominant forearm", "Non-dominant forearm", "Dominant quadriceps", "Non-dominant quadriceps"))
) |>
arrange(site) |> flextable() |>
set_header_labels(
site = "Site", pre_cont = "Pre", post_cont = "Post" , abc_temp_cont = "Diff~post-pre~", pre_exp = "Pre", post_exp = "Post" , abc_temp_exp = "Diff~post-pre~"
) |>
add_header_row(values = c(
"Site", "Control", "Exercise"
), colwidths = c(1,3,3)) |>
colformat_md(part = "all", j = 1:7) |>
merge_v(part = "header") |>
bold(part = "header") |>
align(part = "header", align = "center") |>
autofit() |>
set_table_properties(layout = "autofit") |>
add_footer_lines(values = c("Data from 38 participants. Units: °C"))
save_as_docx(t_summary_temp, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_summary_temp.docx")
t_summary_tempSite | Control | Exercise | ||||
|---|---|---|---|---|---|---|
Pre | Post | Diffpost-pre | Pre | Post | Diffpost-pre | |
Dominant forearm | 31.92 ± 0.96 | 31.42 ± 1.18 | -0.50 ± 0.67 | 31.47 ± 0.90 | 31.21 ± 0.84 | -0.26 ± 0.78 |
Non-dominant forearm | 31.71 ± 0.92 | 31.92 ± 1.04 | 0.21 ± 0.71 | 31.28 ± 0.89 | 32.14 ± 0.80 | 0.86 ± 0.87 |
Dominant quadriceps | 31.30 ± 0.79 | 31.16 ± 1.01 | -0.14 ± 0.62 | 31.03 ± 0.80 | 32.04 ± 0.98 | 1.01 ± 1.01 |
Non-dominant quadriceps | 30.61 ± 1.00 | 30.35 ± 1.12 | -0.26 ± 0.57 | 30.21 ± 0.77 | 31.62 ± 0.92 | 1.41 ± 0.97 |
Data from 38 participants. Units: °C | ||||||
Then we estimate the impact of the factors ‘site’ and ‘condition’ on the post - pre difference in skin temperature while controlling for pre condition values
mtemp <- lmer(abc_temp ~ pre + site*condition + (1|ID), data = delta_temp)
check_mod_s(mtemp, delta_temp)
cooks_outliers(mtemp, delta_temp)No outlier detectedRandom effects slightly deviates from normality. Residuals are homoskedastic and slightly deviates from normality (negative skew) in factors Nondominantforearm:Exercise. Cooks distance < 1 suggest little influence of deviated data point.
summary(mtemp, ddf = "Kenward-Roger")Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: abc_temp ~ pre + site * condition + (1 | ID)
Data: delta_temp
REML criterion at convergence: 682.2
Scaled residuals:
Min 1Q Median 3Q Max
-3.4274 -0.5419 0.0332 0.6751 3.1383
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.1413 0.3758
Residual 0.4320 0.6572
Number of obs: 304, groups: ID, 38
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 10.39375 1.75457 259.02666 5.924 9.97e-09 ***
pre -0.32391 0.05620 259.68767 -5.763 2.33e-08 ***
site1 -0.50524 0.07117 276.00360 -7.099 1.06e-11 ***
site2 0.34138 0.06750 265.81476 5.058 7.92e-07 ***
site3 0.13552 0.06531 258.13435 2.075 0.0390 *
condition1 -0.40079 0.03925 267.35423 -10.211 < 2e-16 ***
site1:condition1 0.35591 0.06531 258.14507 5.449 1.18e-07 ***
site2:condition1 0.14679 0.06530 258.09863 2.248 0.0254 *
site3:condition1 -0.13470 0.06537 258.37352 -2.061 0.0403 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) pre site1 site2 site3 cndtn1 st1:c1 st2:c1
pre -0.999
site1 0.394 -0.394
site2 0.251 -0.251 -0.198
site3 -0.024 0.024 -0.316 -0.329
condition1 0.275 -0.276 0.109 0.069 -0.007
sit1:cndtn1 0.026 -0.026 0.010 0.006 -0.001 0.007
sit2:cndtn1 0.017 -0.017 0.007 0.004 0.000 0.005 -0.333
sit3:cndtn1 -0.050 0.050 -0.020 -0.012 0.001 -0.014 -0.334 -0.334confint(mtemp) 2.5 % 97.5 %
.sig01 0.265004428 0.506409340
.sigma 0.596151233 0.706889857
(Intercept) 6.999364442 13.809437086
pre -0.433321355 -0.215178935
site1 -0.643438161 -0.367564603
site2 0.210367408 0.472007406
site3 0.008955115 0.262124279
condition1 -0.476982926 -0.324834633
site1:condition1 0.229297704 0.482477760
site2:condition1 0.020210461 0.273343357
site3:condition1 -0.261368689 -0.007956001# model anova table, reported
an_temp <- anova(mtemp, ddf = "Kenward-Roger") |>
as_tibble() |>
mutate(om2 = as_tibble(F_to_omega2(`F value`, NumDF, DenDF))[1] |> pull()) |>
mutate(interpretation = as_tibble(interpret_omega_squared(om2))[1] |> pull() )
an_temp$effects <- rownames(anova(mtemp, ddf = "Kenward-Roger"))
t_mtemp <- an_temp |>
rowwise() |>
mutate(`Pr(>F)` = pvalue_format_table(`Pr(>F)`)) |>
select(effects, everything()) |>
flextable() |>
set_header_labels(effects = "Effect", `Sum Sq` = "SS", `Mean Sq` = "MS", NumDF ="Df num", DenDF = "Df den", `F value`= "F value", `Pr(>F)` = "p", om2 = "\u03C9\u00B2\u209A", interpretation = "Interpretation") |>
add_header_row(values = c( "Effect","SS", "MS", "Df num", "Df den", "F value", "p", "Effect size" ), colwidths = c(1,1,1,1,1,1,1,2)
) |>
colformat_double(digits = 3) |>
merge_v(j = c(1,2,3), part = "body") |>
merge_v(part = "header") |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
autofit() |>
add_footer_lines(values = c("Data from 38 participants. \n Abreviations: SS, Sum of squares; MS, Mean squares; Df, Degrees of freedom"))
t_mtempEffect | SS | MS | Df num | Df den | F value | p | Effect size | |
|---|---|---|---|---|---|---|---|---|
ω²ₚ | Interpretation | |||||||
pre | 14.348 | 14.348 | 1 | 259.688 | 33.215 | < .001 | 0.110 | medium |
site | 28.854 | 9.618 | 3 | 270.169 | 22.264 | < .001 | 0.189 | large |
condition | 45.042 | 45.042 | 1 | 267.354 | 104.272 | < .001 | 0.277 | large |
site:condition | 22.920 | 7.640 | 3 | 258.170 | 17.687 | < .001 | 0.160 | large |
Data from 38 participants. | ||||||||
#print(ts6_alt, preview = "docx")
save_as_docx(t_mtemp, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_mtemp.docx")We compute contrasts of marginal means and cohen’s dz based on t value and ddf.
EMM <- emmeans(mtemp, ~site*condition)
pairs.cond <- pairs(EMM, simple = "condition", adjust = "bonferroni") |> as_tibble()
pairs.cond.ci <- pairs(EMM, simple = "condition", adjust = "bonferroni") |>
confint() |>
as_tibble()
pairs.cond <- full_join(pairs.cond, pairs.cond.ci)
pairs.s <- pairs(EMM, simple = "site", adjust = "bonferroni") |> as_tibble()
pairs.s.ci <- pairs(EMM, simple = "site", adjust = "bonferroni") |> confint() |> as_tibble()
pairs.s <- full_join(pairs.s, pairs.s.ci)
pairs_comp <- full_join(
pairs.cond, pairs.s
) |>
rowwise() |>
mutate(d = as_tibble(t_to_d(t.ratio, df, paired = T)[1]) |> pull(),
interpretation = as_tibble(interpret_cohens_d(t_to_d(t.ratio, df, paired = T)))[5] |> pull())
t_cont_temp <- pairs_comp %>%
select(site,condition, contrast, estimate, lower.CL, upper.CL, SE, df, t.ratio, p.value, d, interpretation) %>%
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
flextable() %>%
set_header_labels(site = "Site", condition = "Condition", contrast = "Contrast", estimate = "Difference", lower.CL = "low", upper.CL = "high", SE = "Standard error", df = "Df", t.ratio = "t ratio",p.value = "p", d = "d~z~", interpretation = "interpretation" ) %>% add_header_row(values = c("Site","Condition", "Contrast","Difference","CI 95 %", "Standard error","Df","t ratio", "p", "Effect size" ), colwidths = c(1,1,1,1,2,1,1,1,1,2)) |>
bold(part = "header") %>%
merge_v(part = "header") |>
autofit() %>%
colformat_double( digits = 3) %>%
colformat_md(part = "header", j = 11) |>
fix_border_issues() %>%
set_table_properties(layout = "autofit") |>
merge_v(j = c(1,2,3)) |>
hline(i = c(4), border = dash_border) |>
valign(j = c(1,2,3), valign = "top") |>
bg(i = ~ p.value < 0.05 , j = 10, bg = "yellow" ) |>
add_footer_lines(values = c("Data from 38 participants. \n Units: °C. \n Abreviations: CI 95%, 95% confidence interval; Df, Degrees of freedom; dz , Cohen's d for paired samples."))
t_cont_tempSite | Condition | Contrast | Difference | CI 95 % | Standard error | Df | t ratio | p | Effect size | ||
|---|---|---|---|---|---|---|---|---|---|---|---|
low | high | dz | interpretation | ||||||||
Dominant forearm | Control - Exercise | -0.090 | -0.391 | 0.211 | 0.153 | 261.428 | -0.587 | 0.558 | -0.036 | very small | |
Non-dominant forearm | -0.508 | -0.809 | -0.207 | 0.153 | 261.150 | -3.327 | 0.00101 | -0.206 | small | ||
Dominant quadriceps | -1.071 | -1.369 | -0.773 | 0.152 | 259.326 | -7.067 | < .001 | -0.439 | small | ||
Non-dominant quadriceps | -1.538 | -1.838 | -1.237 | 0.152 | 260.796 | -10.084 | < .001 | -0.624 | medium | ||
Control | Dominant forearm - (Non-dominant forearm) | -0.637 | -1.040 | -0.235 | 0.151 | 258.811 | -4.215 | < .001 | -0.262 | small | |
Dominant forearm - Dominant quadriceps | -0.150 | -0.562 | 0.261 | 0.155 | 264.244 | -0.970 | 1 | -0.060 | very small | ||
Dominant forearm - (Non-dominant quadriceps) | 0.190 | -0.255 | 0.636 | 0.168 | 279.217 | 1.135 | 1 | 0.068 | very small | ||
(Non-dominant forearm) - Dominant quadriceps | 0.487 | 0.082 | 0.893 | 0.153 | 260.900 | 3.195 | 0.00943 | 0.198 | very small | ||
(Non-dominant forearm) - (Non-dominant quadriceps) | 0.828 | 0.395 | 1.261 | 0.163 | 274.490 | 5.082 | < .001 | 0.307 | small | ||
Dominant quadriceps - (Non-dominant quadriceps) | 0.340 | -0.073 | 0.754 | 0.156 | 265.468 | 2.188 | 0.177 | 0.134 | very small | ||
Exercise | Dominant forearm - (Non-dominant forearm) | -1.056 | -1.458 | -0.654 | 0.151 | 258.677 | -6.984 | < .001 | -0.434 | small | |
Dominant forearm - Dominant quadriceps | -1.131 | -1.538 | -0.725 | 0.153 | 261.334 | -7.403 | < .001 | -0.458 | small | ||
Dominant forearm - (Non-dominant quadriceps) | -1.257 | -1.700 | -0.815 | 0.167 | 278.198 | -7.548 | < .001 | -0.453 | small | ||
(Non-dominant forearm) - Dominant quadriceps | -0.076 | -0.478 | 0.327 | 0.151 | 259.163 | -0.499 | 1 | -0.031 | very small | ||
(Non-dominant forearm) - (Non-dominant quadriceps) | -0.202 | -0.633 | 0.230 | 0.162 | 273.886 | -1.243 | 1 | -0.075 | very small | ||
Dominant quadriceps - (Non-dominant quadriceps) | -0.126 | -0.545 | 0.293 | 0.158 | 268.179 | -0.800 | 1 | -0.049 | very small | ||
Data from 38 participants. | |||||||||||
save_as_docx(t_cont_temp, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_cont_temp.docx")Data is visualised:
label <- anova(mtemp, ddf = "Kenward-Roger") |>
tidy() |>
filter(term == "site:condition") |>
mutate(p.value = pvalue_format_plot(p.value)) |>
mutate(interaction = glue("Site x Condition: {p.value}")) |>
pull(interaction)
grob <- grobTree(textGrob(label, x = 0.1, y = 0.95, hjust=0,
gp=gpar(col="black", fontsize=9, fontfamily = "Roboto")))
p_delta_temp_points <- delta_temp %>%
mutate(unique_ID = paste(ID, condition)) %>%
mutate(unique_ID = factor(unique_ID)) %>%
group_by(unique_ID) %>%
arrange(unique_ID) %>%
ggplot() +
geom_point(aes(x = site, y = abc_temp, colour = condition ), position = position_jitter(width = 0.05, seed = 123)) +
geom_line(aes(y = abc_temp, x = site, group = factor(unique_ID), colour = condition), alpha = 0.3, position = position_jitter(width = 0.05, seed = 123)) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
scale_color_manual(values = c( "#00468BFF","#ED0000FF")) +
labs(title = "Skin temperature: absolute change",
subtitle = "Individual values",
y = "Tskin<sub>abc</sub> (°C)",
x = "",
colour = "Condition"
) +
my_theme() +
theme(legend.position = "none",
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()
) +
guides(colour = guide_legend(override.aes = list(size = 4)))
p <- ggplot_build(p_delta_temp_points)
p$data[[1]] <- p$data[[1]] |>
group_by(group) |>
mutate(x2 = if_else(colour == "#00468BFF", x-0.09, x+0.09)) |>
mutate(x = x2) |>
select(-x2) |> ungroup()
p$data[[2]] <- p$data[[2]] |>
group_by(group) |>
mutate(x2 = if_else(colour == "#00468BFF", x-0.09, x+0.09)) |>
mutate(x = x2 ) |>
select(-x2) |> ungroup()
q <- ggplot_gtable(p)
p_delta_temp_points <- as.ggplot(q)
p_delta_temp_mean <- delta_temp_sum |>
mutate(site = factor(site, levels = c(
"Dominant forearm",
"Non-dominant forearm",
"Dominant quadriceps",
"Non-dominant quadriceps"
), labels = c(
"<span style='color:black'>Dominant forearm</span>",
"<span style='color:#925E9FFF'>Non-dominant forearm</span>",
"<span style='color:#42B540FF'>Dominant quadriceps</span>",
"<span style='color:black'>Non-dominant quadriceps</span>"
))) |>
ggplot(aes(x = site, y = m)) +
geom_errorbar(aes( ymin = cil, ymax = cih, group = condition, colour = condition), width = 0.1, position = position_dodge(width=0.2), alpha = 1) +
geom_point(aes( colour = condition, group = condition), position=position_dodge(width=0.2), size = 4) +
geom_line(aes(x = site, y = m, colour = condition, group = condition), position=position_dodge(width=0.2)) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
labs(
subtitle = "Condition means",
y = "Mean Tskin<sub>abc</sub> (°C)",
x = "Site",
colour = "Condition"
) +
my_theme() +
theme(legend.position = "none") +
annotation_custom(grob) +
showSignificance( 2.3, c(0.2, 1), "*", width = -0.05) +
showSignificance( 3.3, c(-0.15, 1.3), "*", width = -0.05) +
showSignificance( 4.3, c(-0.2, 1.6), "*", width = -0.05) +
showSignificance( c(0.9,2), 1.4,"*", width = -0.05, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( c(0.9,3), 1.8,"*", width = -0.05, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF"))+
showSignificance( c(0.9,4), 2.2,"*", width = -0.05, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( c(1,2.1), -1,"*", width = +0.05, segmentParams = list(col = "#00468BFF"), textParams = list(col = "#00468BFF")) +
showSignificance( c(2,3.1), -1.4,"*", width = +0.05, segmentParams = list(col = "#00468BFF"), textParams = list(col = "#00468BFF"))+
showSignificance( c(2,4.1), -1.8,"*", width = +0.05, segmentParams = list(col = "#00468BFF"), textParams = list(col = "#00468BFF")) +
coord_cartesian(ylim = c(-3,4)) +
theme(axis.text.x = element_markdown())
p_delta_temp_mean <- as.ggplot(p_delta_temp_mean)
p_temp <- p_delta_temp_points/p_delta_temp_mean + plot_layout(guides = "collect")
#plot_annotation(tag_levels = 'A')
#ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/ptemp.png", plot = p_temp, height = 9, width = 7.5,dpi = 600)
#p_temp
mean_temp_fun <- function(s){
pos_temp <- position_dodge(width=0.2)
p <- d_temp_sum |>
filter(site == s) |>
ggplot(aes(x = timing, y = m, col = condition, group = condition)) +
geom_point(size = 4, position = pos_temp) +
geom_line(position = pos_temp) +
geom_errorbar(aes(ymin = cil, ymax = cih), position = pos_temp, width = 0.1) +
labs(subtitle = s, x = "", y = "Mean Tskin (°C)", colour = "") +
scale_color_manual(values = c("#00468BFF", "#ED0000FF")) +
ylim(c(30, 34)) +
theme(legend.position = "none") +
my_theme()
return(p)
}
l_temp_plots <- lapply(levels(d_temp_sum$site),mean_temp_fun )
pt1 <- (l_temp_plots[[1]]+ scale_x_discrete(breaks = "") + labs(title = "Skin temperature: pre and post") + theme(
legend.position = c(0.3, 0.8)))
pt2 <- (l_temp_plots[[2]]+ scale_x_discrete(breaks = "") + theme(plot.subtitle = element_text(colour = "#925E9FFF")))
pt3 <- (l_temp_plots[[3]]+ scale_x_discrete(breaks = "")) + theme(plot.subtitle = element_text(colour = "#42B540FF"))
pt4 <- (l_temp_plots[[4]] +labs(x = "Timing"))
p_mean_temp <-( pt1) /( pt2)/ (pt3)/ (pt4) + plot_layout(axis_titles = "collect")
p_mean_temp <- as.ggplot(p_mean_temp)
# fix plot annotation
ptemp_all <- (p_mean_temp) + p_temp + plot_annotation(tag_levels = 'A' ) + plot_layout(width = c(1,2), ncol = 2, nrow = 1)& theme(
plot.tag = element_text(size = 11, family = "Roboto", face = "bold")
)
ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/fig2_temp.png", plot = ptemp_all, height = 11, width = 9.5,dpi = 600)
ptemp_all
4 EIH models
Manipulation checks suggest adequate operationalization of exercise vs control, no carry-over effect, and reliable sensory testing . Differences in skin temperature changes are observed between conditions and testing sites, stressing the importance of including skin temperature as a covariate in EIH models involving thermal test modalities.
We compute statistical model to investigate the differential impact of exercise across test modalities and sites.
4.1 Model: EIH ~ site x test mod
EIH (i.e., (post - pre)Exercise (post-pre)control) is compared across test modalities and testing sites using a linear mixed model with a random effect on the participant while controlling for baseline values (i.e., (preExercise - precontrol)).
# data is formatted for the model
dm <- data_long %>%
select(-session) |>
filter(test_mod != 'adt') |> # removing ADT data
ungroup() |>
rowwise() |> # cdt and ppt are log transformed
mutate(test_value = if_else(test_mod %in% c("cdt", "ppt"), logt(test_value), test_value)) |>
ungroup() |>
group_by(test_mod, site) |>
mutate(test_value = scale(test_value)) %>% # values of each test mod and site are standardized
group_by(site) |>
mutate(temp = scale(temp)) |>
pivot_wider(names_from = c(condition, timing), values_from = c(temp, test_value)) |>
mutate(
eih = (test_value_exp_post - test_value_exp_pre) - (test_value_cont_post - test_value_cont_pre),
test_value_pre = test_value_exp_pre - test_value_cont_pre,
delta_temp = (temp_exp_post - temp_exp_pre) - (temp_cont_post - temp_cont_pre)
) |>
select(ID, test_mod, site, order, eih, test_value_pre, delta_temp) |>
mutate(test_mod = factor(test_mod, levels = c("ppt", "pp","hpt", "cdt"), labels = c("PPT", "PP", "HPT", "CDT")),
site = factor(site, levels = c("forearm", "quad"), labels = c("Forearm", "Quadriceps"))) # empty mod
#empty_mod <- lmer(eih ~ 1 + (1|ID), data = dm)
#empty.mod(empty_mod)
# model reduced vs with temp
m0 <- lmer(eih ~ test_value_pre + test_mod*site + (1|ID), data = dm ) # reduced
check_mod(m0, dm)
#check_model_ac(m0, dm)
cooks_outliers(m0, dm)No outlier detectedRandom effects are normally distributed. Overall, residuals are homoskedastic and normally distributed. A slight deviation from normality for forearm:pp residuals. No outlier is detected.
We compute a second model including the skin temperature (defined as absolute change of skin temperature at each site across interventions (post - pre)Exercise - (post - pre)control).
m1 <- lmer(eih ~ test_value_pre + delta_temp + test_mod*site + (1|ID), data = dm) # + temp
check_mod(m1, dm)
Assumptions checks are similar in both models. The model are compared using LRT.
print(anova(m1, m0))Data: dm
Models:
m0: eih ~ test_value_pre + test_mod * site + (1 | ID)
m1: eih ~ test_value_pre + delta_temp + test_mod * site + (1 | ID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m0 11 2218.3 2271.3 -1098.2 2196.3
m1 12 2220.0 2277.8 -1098.0 2196.0 0.3154 1 0.5744We estimate the reduced model.
# summary table, not reported
ci_m0 <- rownames_to_column(as.data.frame(confint(m0, oldNames = F)), var = "term") %>%
as_tibble() |>
rename("CI_low" = `2.5 %`) |>
rename("CI_high" = `97.5 %`) |>
mutate(term = str_replace_all(term, "sd_\\(Intercept\\)\\|ID", "sd__(Intercept)")) |>
mutate(term = str_replace(term, "sigma", "sd__Observation"))
df_m0 <- m0 |>
tidy(ddf = "Kenward-Roger") |>
rename("p" = "p.value") |>
rename("t" = "statistic") |>
mutate(effect = str_replace_all(effect, "ran_pars", "random"))
ts6 <- full_join(df_m0, ci_m0) %>% # model summary table
as_tibble() |>
filter(effect != "NA") |>
select(effect, term, estimate, std.error, CI_low, CI_high, t, df, p) |>
flextable() %>%
bold(part = "header", bold = TRUE) %>%
set_header_labels( effect = "Effect", term = "Term", estimate = "Estimate", std.error = "Standard error", CI_low = "low", CI_high = "high", t = "t value", df = "Df", p = "p") |>
add_header_row(
values = c( "Effect", "Term","Estimate","Standard error", "CI 95%", "t value", "Df", "p"),
colwidths = c(1, 1, 1, 1, 2, 1,1,1)
) %>%
merge_h(part = "header") %>%
merge_v(part = "header") %>%
align(align = "center", part = "header") %>%
merge_v(part = "body") |>
valign(j = 1, valign = "top") |>
colformat_double(digits = 3) %>%
autofit() %>%
fix_border_issues() %>%
add_footer_lines(values = c(
"CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Pinprick rating; PPT, Pressure pain threshold"
)) %>%
set_table_properties(layout = "autofit")
#print(ts6, preview = "docx")
# model anova table, reported
df_an <- anova(m0, ddf = "Kenward-Roger") |>
as_tibble() |>
mutate(om2 = as_tibble(F_to_omega2(`F value`, NumDF, DenDF))[1] |> pull()) |>
mutate(interpretation = as_tibble(interpret_omega_squared(om2))[1] |> pull() )
df_an$effects <- rownames(anova(m0, ddf = "Kenward-Roger"))
ts6_alt <- df_an |>
rowwise() |>
mutate(`Pr(>F)` = pvalue_format_table(`Pr(>F)`)) |>
select(effects, everything()) |>
flextable() |>
set_header_labels(effects = "Effect", `Sum Sq` = "SS", `Mean Sq` = "MS", NumDF ="Df num", DenDF = "Df den", `F value`= "F value", `Pr(>F)` = "p", om2 = "\u03C9\u00B2\u209A", interpretation = "Interpretation") |>
add_header_row(values = c( "Effect","SS", "MS", "Df num", "Df den", "F value", "p", "Effect size" ), colwidths = c(1,1,1,1,1,1,1,2)
) |>
colformat_double(digits = 3) |>
merge_v(j = c(1,2,3), part = "body") |>
merge_v(part = "header") |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
bg(i = ~ `Pr(>F)`<0.05, j = 7, bg = "yellow") |>
autofit() |>
add_footer_lines(values = c("Data from 38 participants. \n Abreviations: SS, Sum of squares; MS, Mean squares; Df, Degrees of freedom."))
ts6Effect | Term | Estimate | Standard error | CI 95% | t value | Df | p | |
|---|---|---|---|---|---|---|---|---|
low | high | |||||||
fixed | (Intercept) | -0.001 | 0.052 | -0.104 | 0.102 | -0.017 | 36.934 | 0.987 |
test_value_pre | -0.831 | 0.035 | -0.900 | -0.761 | -23.412 | 903.000 | 0.000 | |
test_mod1 | 0.124 | 0.045 | 0.035 | 0.212 | 2.733 | 866.027 | 0.006 | |
test_mod2 | 0.015 | 0.045 | -0.073 | 0.104 | 0.338 | 866.271 | 0.735 | |
test_mod3 | -0.053 | 0.046 | -0.142 | 0.036 | -1.162 | 867.291 | 0.245 | |
site1 | 0.004 | 0.026 | -0.047 | 0.055 | 0.158 | 866.153 | 0.875 | |
test_mod1:site1 | -0.157 | 0.045 | -0.246 | -0.069 | -3.479 | 866.185 | 0.001 | |
test_mod2:site1 | -0.012 | 0.045 | -0.100 | 0.076 | -0.271 | 866.005 | 0.787 | |
test_mod3:site1 | 0.035 | 0.045 | -0.053 | 0.124 | 0.781 | 866.031 | 0.435 | |
random | sd__(Intercept) | 0.276 | 0.199 | 0.369 | ||||
sd__Observation | 0.788 | 0.749 | 0.822 | |||||
CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Pinprick rating; PPT, Pressure pain threshold | ||||||||
ts6_altEffect | SS | MS | Df num | Df den | F value | p | Effect size | |
|---|---|---|---|---|---|---|---|---|
ω²ₚ | Interpretation | |||||||
test_value_pre | 340.195 | 340.195 | 1 | 903.000 | 548.125 | < .001 | 0.377 | large |
test_mod | 5.852 | 1.951 | 3 | 866.546 | 3.143 | 0.0246 | 0.007 | very small |
site | 0.015 | 0.015 | 1 | 866.153 | 0.025 | 0.875 | 0.000 | very small |
test_mod:site | 10.045 | 3.348 | 3 | 866.081 | 5.395 | 0.00111 | 0.015 | small |
Data from 38 participants. | ||||||||
#print(ts6_alt, preview = "docx")
save_as_docx(ts6_alt, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_m.all_anova.docx")4.1.1 Contrasts model site:test mod
We compute simple pairwise comparisons of the significant test modality:site interaction with bonferroni correction.
EMM <- emmeans(m0, ~ test_mod*site)
pairs.tm <- pairs(EMM, simple = "test_mod", adjust = "holm") |> as_tibble()
pairs.tm.ci <- pairs(EMM, simple = "test_mod", adjust = "holm") |>
confint() |>
as_tibble()
pairs.tm <- full_join(pairs.tm, pairs.tm.ci)
pairs.s <- pairs(EMM, simple = "site", adjust = "holm") |> as_tibble()
pairs.s.ci <- pairs(EMM, simple = "site", adjust = "holm") |> confint() |> as_tibble()
pairs.s <- full_join(pairs.s, pairs.s.ci)
pairs_comp <- full_join(
pairs.tm, pairs.s
) |>
rowwise() |>
mutate(d = as_tibble(t_to_d(t.ratio, df, paired = T)[1]) |> pull(),
interpretation = as_tibble(interpret_cohens_d(t_to_d(t.ratio, df, paired = T)))[5] |> pull())
ts7 <- pairs_comp %>%
select(site, test_mod, contrast, estimate, lower.CL, upper.CL, SE, df, t.ratio, p.value, d, interpretation) %>%
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
flextable() %>%
set_header_labels(site = "Site", test_mod = "Test modalities", contrast = "Contrast", estimate = "Difference", lower.CL = "low", upper.CL = "high", SE = "Standard error", df = "Df", t.ratio = "t ratio",p.value = "p", d = "d~z~", interpretation = "interpretation" ) %>% add_header_row(values = c("Site","Test modalities", "Contrast","Difference","CI 95 %", "Standard error","Df","t ratio", "p", "Effect size" ), colwidths = c(1,1,1,1,2,1,1,1,1,2)) |>
bold(part = "header") %>%
merge_v(part = "header") |>
autofit() %>%
colformat_double( digits = 3) %>%
colformat_md(part = "header", j = 11) |>
fix_border_issues() %>%
set_table_properties(layout = "autofit") |>
merge_v(j = c(1,3)) |>
valign(j = 1, valign = "top") |>
bg(i = ~p.value < 0.05, j = 10, bg = "yellow") |>
hline(i = c(6, 12), border = dash_border) |>
add_footer_lines(values = c("Data from 38 participants. \n Units: PPT, log; CDT, log; HPT, °C; PP, VAS 100 - 0. \n Abreviations: PPT, Pressure pain threshold; CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Mechanical pinprick rating; CI 95%, 95% confidence interval; Df, Degrees of freedom; dz , Cohen's d for paired samples."))
ts7Site | Test modalities | Contrast | Difference | CI 95 % | Standard error | Df | t ratio | p | Effect size | ||
|---|---|---|---|---|---|---|---|---|---|---|---|
low | high | dz | interpretation | ||||||||
Forearm | PPT - PP | -0.037 | -0.313 | 0.239 | 0.104 | 866.020 | -0.354 | 1 | -0.012 | very small | |
PPT - HPT | -0.016 | -0.293 | 0.261 | 0.105 | 866.472 | -0.154 | 1 | -0.005 | very small | ||
PPT - CDT | -0.082 | -0.358 | 0.194 | 0.104 | 866.009 | -0.789 | 1 | -0.027 | very small | ||
PP - HPT | 0.021 | -0.256 | 0.298 | 0.105 | 866.671 | 0.198 | 1 | 0.007 | very small | ||
PP - CDT | -0.045 | -0.321 | 0.231 | 0.104 | 866.004 | -0.435 | 1 | -0.015 | very small | ||
HPT - CDT | -0.066 | -0.343 | 0.211 | 0.105 | 866.598 | -0.632 | 1 | -0.021 | very small | ||
Quadriceps | PPT - PP | 0.253 | -0.023 | 0.530 | 0.105 | 866.207 | 2.425 | 0.0621 | 0.082 | very small | |
PPT - HPT | 0.369 | 0.093 | 0.645 | 0.104 | 866.030 | 3.538 | 0.00213 | 0.120 | very small | ||
PPT - CDT | 0.501 | 0.224 | 0.778 | 0.105 | 866.535 | 4.783 | < .001 | 0.162 | very small | ||
PP - HPT | 0.116 | -0.161 | 0.393 | 0.105 | 866.387 | 1.107 | 0.419 | 0.038 | very small | ||
PP - CDT | 0.248 | -0.028 | 0.524 | 0.104 | 866.080 | 2.372 | 0.0621 | 0.081 | very small | ||
HPT - CDT | 0.132 | -0.146 | 0.409 | 0.105 | 866.804 | 1.256 | 0.419 | 0.043 | very small | ||
PPT | Forearm - Quadriceps | -0.307 | -0.511 | -0.102 | 0.104 | 866.033 | -2.937 | 0.0034 | -0.100 | very small | |
PP | -0.016 | -0.221 | 0.189 | 0.104 | 866.022 | -0.155 | 0.877 | -0.005 | very small | ||
HPT | 0.079 | -0.126 | 0.284 | 0.104 | 866.119 | 0.755 | 0.451 | 0.026 | very small | ||
CDT | 0.277 | 0.072 | 0.482 | 0.105 | 866.221 | 2.649 | 0.00822 | 0.090 | very small | ||
Data from 38 participants. | |||||||||||
#print(ts7, preview = "docx")
save_as_docx(ts7, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_m.all_contrasts.docx")4.1.2 Plot model site: test mod
Data are plotted and results from the model are displayed.
txt_gain <- text_grob(label = "Increased sensitivity", size = 9,rot = 90, color = "grey")
txt_loss <- text_grob(label = "Reduced sensitivity", size = 9,rot = 90, color = "grey")
df_label <- anova(m0, ddf = "Kenward-Roger") |> # replaced m0
tidy() |>
rowwise() |>
mutate(p.value = pvalue_format_plot(p.value))
label <- df_label |>
filter(term == "test_mod:site") |>
mutate(interaction = glue("Test modality x Site: {p.value}")) |>
pull(interaction)
grob <- grobTree(textGrob(label, x = 0.1, y = 0.85, hjust=0,
gp=gpar(col="black", fontsize=9, fontfamily = "Roboto")))
data.mod.summary <- dm %>%
group_by(ID, site, test_mod) %>%
mutate(eih = mean(eih)) %>%
slice(1) %>%
ungroup() |>
group_by(test_mod, site) %>%
summarize(m = mean(eih),
sd = sd(eih),
cil = ci_low(eih),
cih = ci_high(eih)
) %>%
mutate(unique_site = paste(test_mod, site)) %>% arrange(unique_site)
pos3 <- position_jitter(width = 0.05, seed = 123)
p.mod.all.points <- dm %>%
group_by(ID, site, test_mod) %>%
mutate(eih = mean(eih)) %>%
slice(1) %>%
ungroup() |>
mutate(unique_ID = paste(ID, site)) %>%
mutate(unique_ID = factor(unique_ID)) %>%
ggplot() +
geom_point(aes(x = test_mod, y = eih, colour = site, group = factor(unique_ID)),position = pos3) +
geom_line(aes(y = eih, x = test_mod, group = factor(unique_ID), colour = site) ,position = pos3, alpha = 0.3) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
scale_color_manual(values = c("#925E9FFF", "#42B540FF")) +
labs(#title = "LMM predicted values",
#subtitle = "EIH ~ Baseline + Test modality * Site",
y = "Individual mean EIH (standardized units)",
x = "",
colour = "Site"
) +
my_theme() +
annotation_custom(txt_gain, ymin = -1.5, ymax = -2, xmin = 0.3, xmax = 0.7) +
annotation_custom(txt_loss, ymin = 1.5, ymax = 2, xmin = 0.3, xmax = 0.7) +
theme(legend.position = "top",
axis.text.x=element_blank(),
axis.ticks.x=element_blank(),
axis.line.x = element_blank()
) + coord_cartesian(ylim = c(-3,4)) +
guides(colour = guide_legend(override.aes = list(size = 4)))
p <- ggplot_build(p.mod.all.points)
p$data[[1]] <- p$data[[1]] |>
group_by(group) |>
mutate(x2 = if_else(colour == "#925E9FFF", x-0.09, x+0.09)) |>
mutate(x = x2) |>
select(-x2) |> ungroup()
p$data[[2]] <- p$data[[2]] |>
group_by(group) |>
mutate(x2 = if_else(colour == "#925E9FFF", x-0.09, x+0.09)) |>
mutate(x = x2 ) |>
select(-x2) |> ungroup()
q <- ggplot_gtable(p)
p.mod.all.points <- as.ggplot(q)
p.mod.all.means <- ggplot() +
geom_errorbar(data =data.mod.summary, aes(x = test_mod, y = m, ymin = cil, ymax = cih, group = site, colour = site), width = 0.1, position = position_dodge(width=0.2), alpha = 1) +
geom_point(data = data.mod.summary, aes(x = test_mod, y = m, colour = site, group = site), position=position_dodge(width=0.2), size = 4) +
geom_line(data = data.mod.summary ,aes(x = test_mod, y = m, colour = site, group = site), position=position_dodge(width=0.2)) +
geom_hline(yintercept = 0, linetype = "dashed", colour = "black") +
scale_color_manual(values = c("#925E9FFF", "#42B540FF")) +
labs(#title = "LMM predicted values",
#subtitle = "EIH ~ Baseline + Test modality * Site",
y = "Condition mean EIH (standardized units)",
x = "Test modality",
colour = "Site"
) +
my_theme() +
theme(legend.position = "none") +
annotation_custom(txt_gain, ymin = -1, ymax = -0.5, xmin = 0.3, xmax = 0.7) +
annotation_custom(txt_loss, ymin = 0.5, ymax = 1, xmin = 0.3, xmax = 0.7) +
annotation_custom(grob)+
coord_cartesian(ylim = c(-1.5,1.5)) +
showSignificance( 0.8,c(-0.18, 0.47) , "*", width = +0.05) +
showSignificance( 4.2, c(-0.55, 0.23), "*", width = -0.05) +
showSignificance( c(1,3.1), 0.45 + 0.6,"*", width = -0.05,textParams = list(col = "#42B540FF"), segmentParams = list(col = "#42B540FF") ) +
showSignificance( c(1,4.1), 0.45 + 0.3,"*", width = -0.05,textParams = list(col = "#42B540FF"), segmentParams = list(col = "#42B540FF"))
p.mod.all.means <- as.ggplot(p.mod.all.means)
p.mod.all <-p.mod.all.points/p.mod.all.means + plot_annotation(tag_levels = 'A') + plot_layout(guides = "collect") & theme(
legend.position = "top",
plot.tag = element_text(size = 11, family = "Roboto", face = "bold")
)
ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/fig3_modall.png", plot = p.mod.all, height = 9, width = 7.5,dpi = 600)
p.mod.all
4.2 Model: absolute change ~ site x condition
d_abc <- data_long |>
select(-session) |>
filter(test_mod != 'adt') |> # removing ADT data
ungroup() |>
rowwise() |> # cdt and ppt are log transformed
mutate(test_value = if_else(test_mod %in% c("cdt", "ppt"), logt(test_value), test_value)) |>
ungroup() |>
group_by(test_mod, site) |>
mutate(test_value = scale(test_value)) %>% # values of each test mod and site are standardized
group_by(site) |>
mutate(temp = scale(temp)) |>
pivot_wider(names_from = timing, values_from = c(test_value, temp)) |>
mutate(
abc = test_value_post - test_value_pre,
pre = test_value_pre,
abc_temp = temp_post - temp_pre
) |>
select(ID, order, condition, site, test_mod, abc, pre, abc_temp)
#. function to filter for each test modality
d.mod.s <- function(tm){
d_abc |> filter(test_mod == tm) |>
mutate(
site = factor(site, levels = c("forearm", "quad")),
condition = factor(condition, levels = c("cont", "exp"))
)
}4.2.1 PPT
dm.ppt <- d.mod.s("ppt") # data set with absolute change
m.0.ppt <- lmer(abc ~ pre + condition*site + (1|ID), data = dm.ppt, REML = T)
check_mod_s(m.0.ppt, dm.ppt)
cooks_outliers(m.0.ppt, dm.ppt)No outlier detectedRandom effects are roughly normally distributed. Residuals are homoskedastic and normally distributed, except for quadriceps:exp. No outlier is detected (cooks distance < 1).
print(summary(m.0.ppt, ddf = "Kenward-Roger"))Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: abc ~ pre + condition * site + (1 | ID)
Data: dm.ppt
REML criterion at convergence: 862
Scaled residuals:
Min 1Q Median 3Q Max
-2.5547 -0.6524 0.0617 0.6773 3.6335
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.2813 0.5304
Residual 0.3020 0.5496
Number of obs: 456, groups: ID, 38
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.004155 0.089805 36.750582 0.046 0.96335
pre -0.541854 0.042321 435.300161 -12.803 < 2e-16 ***
condition1 -0.067884 0.025753 414.169507 -2.636 0.00871 **
site1 -0.005652 0.025736 414.021957 -0.220 0.82628
condition1:site1 0.081086 0.025744 414.087845 3.150 0.00175 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) pre cndtn1 site1
pre 0.001
condition1 0.000 -0.037
site1 0.000 -0.006 0.000
condtn1:st1 0.000 0.025 -0.001 0.000an_ppt <- anova(m.0.ppt, ddf = "Kenward-Roger") |> tidy()
an_ppt$test_mod <- rep("PPT", nrow(an_ppt))
an_ppt |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | test_mod |
|---|---|---|---|---|---|---|---|
| pre | 49.508 | 49.508 | 1 | 435.300 | 163.924 | 0.000 | PPT |
| condition | 2.098 | 2.098 | 1 | 414.170 | 6.948 | 0.009 | PPT |
| site | 0.015 | 0.015 | 1 | 414.022 | 0.048 | 0.826 | PPT |
| condition:site | 2.996 | 2.996 | 1 | 414.088 | 9.921 | 0.002 | PPT |
4.2.2 CDT
dm.cdt <- d.mod.s("cdt") # data set with absolute change
m.0.cdt <- lmer(abc ~ pre + condition*site + (1|ID), data = dm.cdt, REML = T)
check_mod_s(m.0.cdt, dm.cdt)
#check_model_ac(m.0.cdt, dm.cdt)
cooks_outliers(m.0.cdt, dm.cdt)No outlier detectedRandom effects are roughly normally distributed. Residuals are homoskedastic and normally distributed—except for forearm:exp and quad:exp where a slight positive skew is observed. All coods distance < 1: no outlier is detected.
We compute a second model including the skin temperature (defined as absolute change of skin temperature at each site across interventions (post - pre)condition ).
m.1.cdt <- lmer(abc ~ pre + abc_temp + condition*site + (1|ID), data = dm.cdt, REML = T)
check_mod_s(m.1.cdt, dm.cdt)
cooks_outliers(m.1.cdt, dm.cdt)No outlier detectedDistributional assumptions are similar to model without skin temperature absolute change (see above).
print(anova(m.0.cdt, m.1.cdt))Data: dm.cdt
Models:
m.0.cdt: abc ~ pre + condition * site + (1 | ID)
m.1.cdt: abc ~ pre + abc_temp + condition * site + (1 | ID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m.0.cdt 7 967.45 996.3 -476.72 953.45
m.1.cdt 8 968.81 1001.8 -476.41 952.81 0.6318 1 0.4267We choose the reduced model.
print(summary(m.0.cdt, ddf = "Kenward-Roger"))Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: abc ~ pre + condition * site + (1 | ID)
Data: dm.cdt
REML criterion at convergence: 976.7
Scaled residuals:
Min 1Q Median 3Q Max
-2.3150 -0.6578 -0.0244 0.5673 4.4653
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.1952 0.4419
Residual 0.4087 0.6393
Number of obs: 456, groups: ID, 38
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.003022 0.077679 36.678942 -0.039 0.96917
pre -0.501956 0.043814 402.921487 -11.456 < 2e-16 ***
condition1 0.059489 0.030014 414.767332 1.982 0.04813 *
site1 0.058621 0.029986 414.499696 1.955 0.05126 .
condition1:site1 -0.082702 0.029991 414.546350 -2.758 0.00608 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) pre cndtn1 site1
pre -0.001
condition1 0.000 0.071
site1 0.000 0.057 0.004
condtn1:st1 0.000 -0.059 -0.004 -0.003an_cdt <- anova(m.0.cdt, ddf = "Kenward-Roger") |> tidy()
an_cdt$test_mod <- rep("CDT", nrow(an_cdt))
an_cdt |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | test_mod |
|---|---|---|---|---|---|---|---|
| pre | 53.638 | 53.638 | 1 | 402.921 | 131.250 | 0.000 | CDT |
| condition | 1.605 | 1.605 | 1 | 414.767 | 3.929 | 0.048 | CDT |
| site | 1.562 | 1.562 | 1 | 414.500 | 3.822 | 0.051 | CDT |
| condition:site | 3.108 | 3.108 | 1 | 414.546 | 7.604 | 0.006 | CDT |
4.2.3 HPT
dm.hpt <- d.mod.s("hpt") # data set with absolute change
m.0.hpt <- lmer(abc ~ pre + condition*site + (1|ID), data = dm.hpt, REML = T)
check_mod_s(m.0.hpt, dm.hpt)
cooks_outliers(m.0.hpt, dm.hpt)No outlier detectedRandom effects are roughly normally distributed. Residuals are homoskedastic and appear negatively skewed for all site:condition levels. No outlier is detected.
We compute a second model including the skin temperature (defined as absolute change of skin temperature at each site across interventions (post - pre)condition )
m.1.hpt <- lmer(abc ~ pre + abc_temp + condition*site + (1|ID), data = dm.hpt, REML = T)
check_mod_s(m.1.hpt, dm.hpt)
cooks_outliers(m.1.hpt, dm.hpt)No outlier detectedDistributional assumptions are similar to model without skin temperature absolute change (see above).
print(anova(m.0.hpt, m.1.hpt))Data: dm.hpt
Models:
m.0.hpt: abc ~ pre + condition * site + (1 | ID)
m.1.hpt: abc ~ pre + abc_temp + condition * site + (1 | ID)
npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
m.0.hpt 7 832.43 861.28 -409.21 818.43
m.1.hpt 8 828.44 861.42 -406.22 812.44 5.9889 1 0.0144 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1We choose the model including skin temperature model.
print(summary(m.1.hpt, ddf = "Kenward-Roger"))Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: abc ~ pre + abc_temp + condition * site + (1 | ID)
Data: dm.hpt
REML criterion at convergence: 842
Scaled residuals:
Min 1Q Median 3Q Max
-4.1844 -0.5642 0.0706 0.6378 2.3225
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.1224 0.3499
Residual 0.3043 0.5516
Number of obs: 456, groups: ID, 38
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.02570 0.06519 42.59164 0.394 0.6954
pre -0.43764 0.03874 310.51138 -11.298 <2e-16 ***
abc_temp -0.09440 0.03777 448.86594 -2.499 0.0128 *
condition1 -0.05559 0.03134 435.89343 -1.774 0.0768 .
site1 0.03592 0.02592 413.71290 1.386 0.1666
condition1:site1 -0.01927 0.02633 416.28224 -0.732 0.4648
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) pre abc_tm cndtn1 site1
pre -0.032
abc_temp -0.290 0.083
condition1 -0.159 -0.073 0.552
site1 0.022 0.023 -0.076 -0.045
condtn1:st1 0.055 -0.059 -0.189 -0.099 0.013an_hpt <- anova(m.1.hpt, ddf = "Kenward-Roger") |> tidy()
an_hpt$test_mod <- rep("HPT", nrow(an_hpt))
an_hpt |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | test_mod |
|---|---|---|---|---|---|---|---|
| pre | 38.842 | 38.842 | 1 | 310.511 | 127.653 | 0.000 | HPT |
| abc_temp | 1.901 | 1.901 | 1 | 448.866 | 6.247 | 0.013 | HPT |
| condition | 0.957 | 0.957 | 1 | 435.893 | 3.146 | 0.077 | HPT |
| site | 0.584 | 0.584 | 1 | 413.713 | 1.920 | 0.167 | HPT |
| condition:site | 0.163 | 0.163 | 1 | 416.282 | 0.535 | 0.465 | HPT |
4.2.4 PP
dm.pp <- d.mod.s("pp") # data set with absolute change
m.0.pp <- lmer(abc ~ pre + condition*site + (1|ID), data = dm.pp, REML = T)
check_mod_s(m.0.pp, dm.pp)
#check_model_ac(m.0.pp, dm.pp)
cooks_outliers(m.0.pp, dm.pp)No outlier detectedRandom effects are roughly normally distributed. Residuals are homoskedastic and follow a normal distribution with fat tails. No outlier is detected (cooks distance < 1).
print(summary(m.0.pp, ddf = "Kenward-Roger"))Linear mixed model fit by REML. t-tests use Kenward-Roger's method [
lmerModLmerTest]
Formula: abc ~ pre + condition * site + (1 | ID)
Data: dm.pp
REML criterion at convergence: 906.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.00847 -0.53594 0.00263 0.47718 3.08691
Random effects:
Groups Name Variance Std.Dev.
ID (Intercept) 0.3869 0.6220
Residual 0.3275 0.5722
Number of obs: 456, groups: ID, 38
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.064036 0.104427 36.758084 -0.613 0.544
pre -0.644480 0.050121 426.162263 -12.859 <2e-16 ***
condition1 -0.001871 0.026837 414.369751 -0.070 0.944
site1 -0.001242 0.026798 414.021745 -0.046 0.963
condition1:site1 0.001728 0.026805 414.087397 0.064 0.949
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) pre cndtn1 site1
pre -0.023
condition1 -0.001 0.053
site1 0.000 -0.002 0.000
condtn1:st1 0.001 -0.023 -0.001 0.000an_pp <- anova(m.0.pp, ddf = "Kenward-Roger") |> tidy()
an_pp$test_mod <- rep("PP", nrow(an_pp))
an_pp |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | test_mod |
|---|---|---|---|---|---|---|---|
| pre | 54.144 | 54.144 | 1 | 426.162 | 165.343 | 0.000 | PP |
| condition | 0.002 | 0.002 | 1 | 414.370 | 0.005 | 0.944 | PP |
| site | 0.001 | 0.001 | 1 | 414.022 | 0.002 | 0.963 | PP |
| condition:site | 0.001 | 0.001 | 1 | 414.087 | 0.004 | 0.949 | PP |
4.2.5 Table models condition:site
F tests and p values on fixed and interactions effects of models are collated in a table.
an_mcs <- rbind(an_ppt, an_cdt, an_hpt, an_pp) |>
select(test_mod, everything()) |>
rowwise() |>
mutate(om2 = as_tibble(F_to_omega2(statistic, NumDF, DenDF))[1] |> pull()) |>
mutate(interpretation = as_tibble(interpret_omega_squared(om2))[1] |> pull() ) |>
mutate(test_mod = factor(test_mod, levels = c("CDT", "HPT", "PP", "PPT")))
t_models_cs <- an_mcs |>
select(test_mod, term, sumsq, meansq, NumDF, DenDF, statistic, p.value, om2, interpretation) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
flextable() |>
set_header_labels( test_mod = "Test modalities", term = "Effect", sumsq = "SS", meansq = "MS", NumDF ="Df num", DenDF = "Df den", statistic= "F value", p.value = "p", om2 = "\u03C9\u00B2\u209A", interpretation = "Interpretation") |>
add_header_row(values = c( "Test modalities", "Effect","SS", "MS", "Df num", "Df den", "F value", "p", "Effect size" ), colwidths = c(1,1,1,1,1,1,1,1,2)
) |> colformat_double(digits = 3) |>
merge_v(j = c(1,2), part = "body") |>
merge_v(part = "header") |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
bg(i = ~p.value < 0.05, j = 8, bg = "yellow") |>
autofit() |>
hline(i = c(4, 8, 13), border = dash_border) |>
add_footer_lines(values = c("Data from 38 participants. \n Units: PPT, log; CDT, log; HPT, °C; PP, VAS 100 - 0; ADT, mA. \n Abreviations: PPT, Pressure pain threshold; CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Mechanical pinprick rating; SS, Sum of squares; MS, Mean squares; Df, Degrees of freedom"))
#print(t_models_pp, preview = "docx")
t_models_csTest modalities | Effect | SS | MS | Df num | Df den | F value | p | Effect size | |
|---|---|---|---|---|---|---|---|---|---|
ω²ₚ | Interpretation | ||||||||
PPT | pre | 49.508 | 49.508 | 1 | 435.300 | 163.924 | < .001 | 0.271 | large |
condition | 2.098 | 2.098 | 1 | 414.170 | 6.948 | 0.00871 | 0.014 | small | |
site | 0.015 | 0.015 | 1 | 414.022 | 0.048 | 0.826 | 0.000 | very small | |
condition:site | 2.996 | 2.996 | 1 | 414.088 | 9.921 | 0.00175 | 0.021 | small | |
CDT | pre | 53.638 | 53.638 | 1 | 402.921 | 131.250 | < .001 | 0.243 | large |
condition | 1.605 | 1.605 | 1 | 414.767 | 3.929 | 0.0481 | 0.007 | very small | |
site | 1.562 | 1.562 | 1 | 414.500 | 3.822 | 0.0513 | 0.007 | very small | |
condition:site | 3.108 | 3.108 | 1 | 414.546 | 7.604 | 0.00608 | 0.016 | small | |
HPT | pre | 38.842 | 38.842 | 1 | 310.511 | 127.653 | < .001 | 0.288 | large |
abc_temp | 1.901 | 1.901 | 1 | 448.866 | 6.247 | 0.0128 | 0.012 | small | |
condition | 0.957 | 0.957 | 1 | 435.893 | 3.146 | 0.0768 | 0.005 | very small | |
site | 0.584 | 0.584 | 1 | 413.713 | 1.920 | 0.167 | 0.002 | very small | |
condition:site | 0.163 | 0.163 | 1 | 416.282 | 0.535 | 0.465 | 0.000 | very small | |
PP | pre | 54.144 | 54.144 | 1 | 426.162 | 165.343 | < .001 | 0.277 | large |
condition | 0.002 | 0.002 | 1 | 414.370 | 0.005 | 0.944 | 0.000 | very small | |
site | 0.001 | 0.001 | 1 | 414.022 | 0.002 | 0.963 | 0.000 | very small | |
condition:site | 0.001 | 0.001 | 1 | 414.087 | 0.004 | 0.949 | 0.000 | very small | |
Data from 38 participants. | |||||||||
save_as_docx(t_models_cs, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_models_cs.docx")4.2.6 Pairwise comparisons models condition:site
Pairswise comparisons of significant main effects or interactions are computed.
fun_emm_cs <- function(mod, tm.name){
emm <- emmeans(mod, ~site*condition)
pairs.s <- pairs(emm, simple = "site", adjust = "holm") |> as_tibble()
pairs.s.ci <- pairs(emm, simple = "site", adjust = "holm") |>
confint() |> as_tibble()
pairs.s.full <- full_join(pairs.s, pairs.s.ci)
pairs.c <- pairs(emm, simple = "condition", adjust = "holm") |> as_tibble()
pairs.c.ci <- pairs(emm, simple = "condition", adjust = "holm") |>
confint() |> as_tibble()
pairs.c.full <- full_join(pairs.c, pairs.c.ci)
pairs <- full_join(pairs.s.full, pairs.c.full) |>
select(condition, site, contrast, estimate, lower.CL, upper.CL, everything())
pairs$var <- rep(tm.name, nrow(pairs))
return(pairs)
} # function for pairwise comparisons of significant interactions
pairs.ppt <- fun_emm_cs(m.0.ppt, "PPT")
pairs.cdt <- fun_emm_cs(m.0.cdt, "CDT")
df_cont_mcs <- rbind(pairs.ppt, pairs.cdt) |>
rename(test_mod = var)
#df_cont_mpp$p.corr <- p.adjust(df_cont_mpp$p.value, method = "BH" )
ts_mcs_contrasts <- df_cont_mcs |>
rowwise() |>
mutate(d = as_tibble(t_to_d(t.ratio, df, paired = T)[1]) |> pull(),
interpretation = as_tibble(interpret_cohens_d(t_to_d(t.ratio, df, paired = T)))[5] |> pull()) |>
select(test_mod, site, condition,contrast, estimate, lower.CL, upper.CL, everything()) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
mutate(test_mod = factor(test_mod, levels = c("CDT", "HPT", "PP", "PPT"))) |>
arrange(test_mod) |>
flextable() |>
set_header_labels(test_mod = "Test modalities", site = "Site", condition = "Condition", contrast = "Contrast", estimate = "Difference", lower.CL = "low", upper.CL = "high", SE = "SE", df = "Df", t.ratio = "t value", p.value = "p", d = "d~z~", interpretation = "Interpretation") |>
add_header_row(values = c( "Test modalities", "Site","Condition","Contrast","Difference","CI", "SE", "Df", "t value", "p", "Effect size" ), colwidths = c(1,1,1,1,1,2,1,1,1,1,2)
) |>
colformat_double(digits = 3) |>
colformat_md(part = "header", j = 12) |>
merge_v(part = "header") |>
merge_v(j = 1) |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
hline(i = 4, border = dash_border) |>
bg(i = ~p.value < 0.05, j = 11, bg = "yellow") |>
autofit() |>
add_footer_lines(values = c("Data from 38 participants. \n Units: PPT, log; CDT, log; HPT, °C. \n Abreviations: PPT, Pressure pain threshold; CDT, Cold detection threshold; HPT, Heat pain threshold; SE, Standard error; CI 95%, 95% confidence interval; Df, Degrees of freedom; dz, Cohen's d for paired samples. "))
#print(ts_mpp_contrasts, preview = "docx")
ts_mcs_contrastsTest modalities | Site | Condition | Contrast | Difference | CI | SE | Df | t value | p | Effect size | ||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
low | high | dz | Interpretation | |||||||||
CDT | cont | forearm - quad | -0.048 | -0.215 | 0.118 | 0.085 | 414.031 | -0.569 | 0.57 | -0.028 | very small | |
exp | forearm - quad | 0.283 | 0.116 | 0.450 | 0.085 | 415.009 | 3.327 | < .001 | 0.163 | very small | ||
forearm | cont - exp | -0.046 | -0.213 | 0.120 | 0.085 | 414.040 | -0.548 | 0.584 | -0.027 | very small | ||
quad | cont - exp | 0.284 | 0.117 | 0.452 | 0.085 | 415.265 | 3.344 | < .001 | 0.164 | very small | ||
PPT | cont | forearm - quad | 0.151 | 0.008 | 0.294 | 0.073 | 414.037 | 2.072 | 0.0388 | 0.102 | very small | |
exp | forearm - quad | -0.173 | -0.317 | -0.030 | 0.073 | 414.073 | -2.383 | 0.0176 | -0.117 | very small | ||
forearm | cont - exp | 0.026 | -0.117 | 0.169 | 0.073 | 414.025 | 0.363 | 0.717 | 0.018 | very small | ||
quad | cont - exp | -0.298 | -0.441 | -0.155 | 0.073 | 414.232 | -4.089 | < .001 | -0.201 | small | ||
Data from 38 participants. | ||||||||||||
save_as_docx(ts_mcs_contrasts, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/ts_mcs_contrasts.docx")4.2.7 Plot models condtion:site
pppt.points <- plot_eih_s(dm.ppt, t = "Pressure pain threshold", u = "Log", mod = m.0.ppt, x = F)[[2]]
pppt.means <- plot_eih_s(dm.ppt, t = "Pressure pain threshold", u = "Log", mod = m.0.ppt, x = F, y_axis = F)[[3]]+
showSignificance( 2.2, c(-0.25, 0.25), "*", width = -0.05) +
showSignificance( c(1.05,2.1), 0.4, "*", width = -0.05, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) + showSignificance( c(0.9,1.95), -0.4, "*", width = 0.05, segmentParams = list(col = "#00468BFF"), textParams = list(col = "#00468BFF")) + theme(legend.position = "none")
pcdt.points <- plot_eih_s(dm.cdt, t = "Cold detection threshold", u = "Log", mod = m.0.cdt, x = T)[[2]]
pcdt.means <- plot_eih_s(dm.cdt, t = "Cold detection threshold", u = "Log", mod = m.0.cdt,x = T, y_axis = F)[[3]] + showSignificance( 2.2, c(-0.4, 0.25), "*", width = -0.05) +
showSignificance( c(1.05,2.1), -0.6, "*", width = 0.05, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF"))+ theme(legend.position = "none")
phpt.points <- plot_eih_s(dm.hpt, t = "Heat pain threshold", u = "°C", mod = m.0.hpt)[[2]]
phpt.means <- plot_eih_s(dm.hpt, t = "Heat pain threshold", u = "°C", mod = m.0.hpt, y_axis = F)[[3]]+ theme(legend.position = "none")
ppp.points <- plot_eih_s(dm.pp, t = "Pinprick rating", u = "VAS 0 - 100", mod = m.0.pp)[[2]]
ppp.means <- plot_eih_s(dm.pp, t = "Pinprick rating", u = "VAS 0 - 100", mod = m.0.pp, y_axis = F)[[3]]+ theme(legend.position = "none")
p_mcs <- (pppt.points + pppt.means) / (ppp.points + ppp.means) / (phpt.points + phpt.means) / (pcdt.points + pcdt.means) +
# plot_layout(ncol = 3, nrow = 4, byrow = T, guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(
# legend.position = "bottom",
plot.tag = element_text(size = 11, family = "Roboto", face = "bold")
)
ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/fig4_mcs.png", plot = p_mcs , dpi = 600, height = 12, width = 10)
p_mcs
4.3 Models value ~ timing * condition | test modality:site
The effects of timing condition on sensory test values are investigated in each test modaltiy:site where a significant of condition was found, and in ADT, using linear mixed model with a random intercept on the participant.
4.3.1 PPT
4.3.1.1 Quadriceps
# input: test modality and site
# output: data set formatted to be used by a lme : test_value ~ condition*timing + (1|ID)
# output data set is in long format and contains all covariables centered and reduced.
d.mod.tms <- function(test_mod, site){
# function that formats data to use in a lme, filtering by site and test modality
# modular output for data accoding to condition or fully prep for lme
tm <- toupper(as.character(test_mod))
s <- switch(site, "forearm" = "Forearm", "quad" = "Quadriceps", "nosite" = "No site")
## data
d.mod <- data_long %>%
select(ID, order, test_mod, site, timing, condition, temp, test_value) |>
mutate(test_mod = factor(test_mod, levels = c("cdt", "hpt", "adt","pp", "ppt"), labels = c("CDT", "HPT", "ADT","PP", "PPT")),
site = factor(site, levels = c("forearm", "quad", "nosite"), labels = c("Forearm", "Quadriceps", "No site")))
if(tm == "ADT"){
d.test <- d.mod |> filter(test_mod == "ADT")
} else {
d.test <- d.mod |>
filter(test_mod ==tm, site == s) |>
ungroup() |>
mutate(temp = scale(temp)) |>
ungroup() |>
mutate(timing = factor(timing, levels = c("pre", "post")),
condition = factor(condition, levels = c("cont", "exp")),
order = factor(order))
}
return(d.test)
}dm.ppt.quad <- d.mod.tms("ppt","quad") |>
rowwise() |>
mutate(test_value = logt(test_value))
#dm.ppt.quad |>
# ggplot(aes(x = timing, y = test_value, col = condition)) +
# geom_boxplot() +
# facet_grid(cols = vars(condition))
m.0.pptquad <- lmer(test_value ~ timing*condition + (1|ID), data = dm.ppt.quad, REML = T)
cooks_outliers( m.0.pptquad, dm.ppt.quad)No outlier detectedcheck_mod_tms(m.0.pptquad, dm.ppt.quad)
Random effects are normally distributed. Residuals are homoskedastics and normally distributed, a slight deviation is observed for post exp (negative skew). No outlier is detected. Cooks distance < 1 suggest little influence of deviated data point. We keep the model.
an_pptq <- anova(m.0.pptquad, ddf = "Kenward-Roger") |> tidy()
an_pptq$site <- rep("Quadriceps", 3)
an_pptq$test_mod <- rep("PPT", 3)
an_pptq |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | site | test_mod |
|---|---|---|---|---|---|---|---|---|
| timing | 0.000 | 0.000 | 1 | 415 | 0.079 | 0.779 | Quadriceps | PPT |
| condition | 0.018 | 0.018 | 1 | 415 | 3.821 | 0.051 | Quadriceps | PPT |
| timing:condition | 0.058 | 0.058 | 1 | 415 | 12.592 | 0.000 | Quadriceps | PPT |
#plot_prepost(dm.ppt.quad , "PPT", "PPT<sub>log</sub> quadriceps ", m.0.pptquad) 4.3.2 CDT
4.3.2.1 Quadriceps
dm.cdt.quad <- d.mod.tms("cdt","quad") |>
rowwise() |>
mutate(test_value = logt(test_value))
m.0.cdtquad <- lmer(test_value ~ timing*condition + (1|ID), data = dm.cdt.quad)
check_mod_tms(m.0.cdtquad, dm.cdt.quad)
An outlier is detected and removed.
dm.cdt.quad.noout <- cooks_outliers( m.0.cdtquad, dm.cdt.quad)1 outliers detected. Observation(s) number: 105m.0.cdtquad.noout <- lmer(test_value ~ timing*condition + (1|ID), data = dm.cdt.quad.noout)
check_mod_tms(m.0.cdtquad.noout, dm.cdt.quad.noout)
Random effects are not fully normal (positvie skew). Residuals are homoskedastics and normal. No further outlier is detected. Data was log transformed according to standard procedure. We decide to keep the model as it is.
an_cdtq <- anova(m.0.cdtquad.noout, ddf = "Kenward-Roger") |> tidy()
an_cdtq$site <- rep("Quadriceps",nrow(an_cdtq))
an_cdtq$test_mod <- rep("CDT", nrow(an_cdtq))
an_cdtq |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | site | test_mod |
|---|---|---|---|---|---|---|---|---|
| timing | 0.050 | 0.050 | 1 | 414.007 | 3.655 | 0.057 | Quadriceps | CDT |
| condition | 0.003 | 0.003 | 1 | 414.007 | 0.244 | 0.621 | Quadriceps | CDT |
| timing:condition | 0.137 | 0.137 | 1 | 414.007 | 10.079 | 0.002 | Quadriceps | CDT |
#plot_prepost(dm.cdt.quad.noout, "CDT", "CDT<sub>log</sub> quadriceps (|\u0394T|)", m.1.cdtquad.noout) 4.3.3 ADT
dm.adt <- d.mod.tms("adt","no site") |>
mutate(condition = factor(condition, levels = c("cont", "exp"))) |>
select(!temp, -site)
m.0.adt <- lmer(test_value ~ timing*condition + (1|ID), data = dm.adt)
cooks_outliers(m.0.adt, dm.adt)No outlier detectedcheck_mod_tms(m.0.adt, na.omit(dm.adt) )
Random effects are not normally distributed. Residuals are homoskedastic and appear normal with fat tails. Extreme values appear dispaly cooks distance < 1 suggesting little impact of models output.
We apply a log transformation.
dm.adt.log <- d.mod.tms("adt","no site") |>
mutate(condition = factor(condition, levels = c("cont", "exp"))) |>
select(!temp, -site) |>
rowwise() |>
mutate(test_value = logt(test_value))
m.0.adt.log <- lmer(test_value ~ timing*condition + (1|ID), data = dm.adt.log)
cooks_outliers(m.0.adt.log, dm.adt.log)No outlier detectedcheck_mod_tms(m.0.adt.log, na.omit(dm.adt.log) )
No impact of log transformation is observed on the distribution of random effects. We decide to estimate the model with orignial units.
an_adt <- anova(m.0.adt, ddf = "Kenward-Roger") |> tidy()
an_adt$site <- rep("No site", nrow(an_adt))
an_adt$test_mod <- rep("ADT", nrow(an_adt))
an_adt |>
kbl(digits = 3) |>
kable_minimal()| term | sumsq | meansq | NumDF | DenDF | statistic | p.value | site | test_mod |
|---|---|---|---|---|---|---|---|---|
| timing | 0 | 0 | 1 | 412.44 | 0.001 | 0.975 | No site | ADT |
| condition | 0 | 0 | 1 | 412.44 | 1.937 | 0.165 | No site | ADT |
| timing:condition | 0 | 0 | 1 | 412.44 | 1.106 | 0.294 | No site | ADT |
#plot_prepost(dm.adt, "ADT", "ADT (mA)", m.0.adt)4.3.4 Table models timing:condition
an_mpp <- rbind(an_pptq, an_cdtq, an_adt) |>
select(test_mod, site, everything()) |>
rowwise() |>
mutate(om2 = as_tibble(F_to_omega2(statistic, NumDF, DenDF))[1] |> pull()) |>
mutate(interpretation = as_tibble(interpret_omega_squared(om2))[1] |> pull() ) |>
mutate(var2 = glue("{test_mod}_{site}"))
t_models_pp <- an_mpp |>
select(test_mod, site, term, sumsq, meansq, NumDF, DenDF, statistic, p.value, om2, interpretation) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value))|>
flextable() |>
set_header_labels( test_mod = "Test modalities", site = "Site", term = "Effect", sumsq = "SS", meansq = "MS", NumDF ="Df num", DenDF = "Df den", statistic= "F value", p.value = "p", om2 = "\u03C9\u00B2\u209A", interpretation = "Interpretation") |>
add_header_row(values = c( "Test modalities", "Site","Effect","SS", "MS", "Df num", "Df den", "F value", "p", "Effect size" ), colwidths = c(1,1,1,1,1,1,1,1,1,2)
) |> colformat_double(digits = 3) |>
hline(i = c(3, 6), border = dash_border) |>
merge_v(j = c(1,2,3), part = "body") |>
merge_v(part = "header") |>
valign(part = "body", valign = "top") |>
bold(part = "header") |>
autofit() |>
bg(i = ~p.value < 0.05, j = 9, bg = "yellow") |>
add_footer_lines(values = c("Data from 38 participants. \n Units: PPT, log; CDT, log; ADT, mA. \n Abreviations: PPT, Pressure pain threshold; CDT, Cold detection threshold; ADT, Auditory detection threshold; SS, Sum of squares; MS, Mean squares; Df, Degrees of freedom"))
#print(t_models_pp, preview = "docx")
t_models_pp Test modalities | Site | Effect | SS | MS | Df num | Df den | F value | p | Effect size | |
|---|---|---|---|---|---|---|---|---|---|---|
ω²ₚ | Interpretation | |||||||||
PPT | Quadriceps | timing | 0.000 | 0.000 | 1 | 415.000 | 0.079 | 0.779 | 0.000 | very small |
condition | 0.018 | 0.018 | 1 | 415.000 | 3.821 | 0.0513 | 0.007 | very small | ||
timing:condition | 0.058 | 0.058 | 1 | 415.000 | 12.592 | < .001 | 0.027 | small | ||
CDT | timing | 0.050 | 0.050 | 1 | 414.007 | 3.655 | 0.0566 | 0.006 | very small | |
condition | 0.003 | 0.003 | 1 | 414.007 | 0.244 | 0.621 | 0.000 | very small | ||
timing:condition | 0.137 | 0.137 | 1 | 414.007 | 10.079 | 0.00161 | 0.021 | small | ||
ADT | No site | timing | 0.000 | 0.000 | 1 | 412.440 | 0.001 | 0.975 | 0.000 | very small |
condition | 0.000 | 0.000 | 1 | 412.440 | 1.937 | 0.165 | 0.002 | very small | ||
timing:condition | 0.000 | 0.000 | 1 | 412.440 | 1.106 | 0.294 | 0.000 | very small | ||
Data from 38 participants. | ||||||||||
save_as_docx(t_models_pp, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t_models_pp.docx")4.3.5 Pairwise comparisons timing:condition
Simple pairwise comparisons of significant main effects and interaction are collated in a table.
# PPT Q
EMM.pptq <- emmeans(m.0.pptquad, ~ timing*condition)
pairs.pptqt <- pairs(EMM.pptq , simple = "timing", adjust = "bonferroni") |> as_tibble()
pairs.pptqt.ci <- pairs(EMM.pptq , simple = "timing", adjust = "bonferroni") |>
confint() |> as_tibble()
pairs.pptqt <- full_join(pairs.pptqt,pairs.pptqt.ci )
pairs.pptqc <- pairs(EMM.pptq, simple = "condition", adjust = "bonferroni") |> as_tibble() # compare doses for each treat
pairs.pptqc.ci <- pairs(EMM.pptq, simple = "condition", adjust = "bonferroni") |>
confint() |>
as_tibble()
pairs.pptqc <- full_join(pairs.pptqc,pairs.pptqc.ci )
pairs.pptq <- full_join(pairs.pptqt, pairs.pptqc)
pairs.pptq$site <- rep("Quadriceps", nrow(pairs.pptq))
pairs.pptq$test_mod <- rep("PPT", nrow(pairs.pptq))
# CDT Q
EMM.cdtq <- emmeans(m.0.cdtquad.noout, ~ timing*condition)
pairs.cdtqt <- pairs(EMM.cdtq , simple = "timing", adjust = "bonferroni") |> as_tibble()
pairs.cdtqt.ci <- pairs(EMM.cdtq , simple = "timing", adjust = "bonferroni") |>
confint() |> as_tibble()
pairs.cdtqt <- full_join(pairs.cdtqt,pairs.cdtqt.ci )
pairs.cdtqc <- pairs(EMM.cdtq, simple = "condition", adjust = "bonferroni") |> as_tibble() # compare doses for each treat
pairs.cdtqc.ci <- pairs(EMM.cdtq, simple = "condition", adjust = "bonferroni") |>
confint() |>
as_tibble()
pairs.cdtqc <- full_join(pairs.cdtqc,pairs.cdtqc.ci )
pairs.cdtq <- full_join(pairs.cdtqt, pairs.cdtqc)
pairs.cdtq$site <- rep("Quadriceps", nrow(pairs.cdtq))
pairs.cdtq$test_mod <- rep("CDT", nrow(pairs.cdtq))
df_cont_mpp <- rbind(pairs.pptq,pairs.cdtq)
ts_mpp_contrasts <- df_cont_mpp |>
rowwise() |>
mutate(d = as_tibble(t_to_d(t.ratio, df, paired = T)[1]) |> pull(),
interpretation = as_tibble(interpret_cohens_d(t_to_d(t.ratio, df, paired = T)))[5] |> pull()) |>
select(test_mod, site, timing, condition,contrast, estimate, lower.CL, upper.CL, everything()) |>
rowwise() |>
mutate(p.value = pvalue_format_table(p.value)) |>
flextable() |>
set_header_labels(test_mod = "Test modalities", site = "Site", timing = "Timing", condition = "Condition", contrast = "Contrast", estimate = "Difference", lower.CL = "low", upper.CL = "high", SE = "SE", df = "Df", t.ratio = "t value", p.value = "p", d = "d~z~", interpretation = "Interpretation") |>
add_header_row(values = c( "Test modalities", "Site","Main effects","Contrast","Difference","CI", "SE", "Df", "t value", "p", "Effect size" ), colwidths = c(1,1,2,1,1,2,1,1,1,1,2)
) |>
colformat_double(digits = 3) |>
colformat_md(part = "header", j = 13) |>
merge_v(part = "header") |>
merge_v(j = c(1,2)) |>
valign(part = "body", valign = "top") |>
hline(i = 4, border = dash_border) |>
bold(part = "header") |>
bg(i = ~p.value < 0.05, j = 12, bg = "yellow") |>
autofit() |>
add_footer_lines(values = c("Data from 38 participants. \n Units: PPT, log; CDT, log \n Abreviations: PPT, Pressure pain threshold; CDT, Cold detection threshold; SE, Standard error; CI 95%, 95% confidence interval; Df, Degrees of freedom; dz, Cohen's d for paired samples. "))
#print(ts_mpp_contrasts, preview = "docx")
ts_mpp_contrastsTest modalities | Site | Main effects | Contrast | Difference | CI | SE | Df | t value | p | Effect size | |||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Timing | Condition | low | high | dz | Interpretation | ||||||||
PPT | Quadriceps | cont | pre - post | 0.021 | 0.003 | 0.039 | 0.009 | 415.000 | 2.310 | 0.0214 | 0.113 | very small | |
exp | pre - post | -0.024 | -0.042 | -0.007 | 0.009 | 415.000 | -2.708 | 0.00704 | -0.133 | very small | |||
pre | cont - exp | 0.010 | -0.008 | 0.028 | 0.009 | 415.000 | 1.127 | 0.26 | 0.055 | very small | |||
post | cont - exp | -0.035 | -0.053 | -0.017 | 0.009 | 415.000 | -3.891 | < .001 | -0.191 | very small | |||
CDT | cont | pre - post | -0.014 | -0.044 | 0.017 | 0.015 | 414.014 | -0.892 | 0.373 | -0.044 | very small | ||
exp | pre - post | 0.056 | 0.025 | 0.086 | 0.015 | 414.000 | 3.601 | < .001 | 0.177 | very small | |||
pre | cont - exp | -0.029 | -0.060 | 0.001 | 0.015 | 414.014 | -1.893 | 0.059 | -0.093 | very small | |||
post | cont - exp | 0.040 | 0.010 | 0.070 | 0.015 | 414.000 | 2.598 | 0.00972 | 0.128 | very small | |||
Data from 38 participants. | |||||||||||||
save_as_docx(ts_mpp_contrasts, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/ts_mpp_contrasts.docx")4.3.6 Plots models timing:condition
# PPT Q
d.sum.pptq <- dm.ppt.quad |>
group_by(timing, condition, ID) |>
mutate(tm = mean(test_value, na.rm = T)) |>
slice(1) |>
ungroup() |>
group_by(timing, condition) |>
summarise(m = mean(tm, na.rm = T),
cil = ci_low(tm),
cih = ci_high(tm))
ppptq <- plot_prepost_an(dm.ppt.quad, t = "Pressure pain threshold", u = "Log", df_an = an_mpp, v = "PPT_Quadriceps", y = "test_value", x = F) + showSignificance( c(1.05,2.1), max(d.sum.pptq$cih) +0.001*max(d.sum.pptq$cih) + 0.01, "*", width = -0.005, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( c(0.9,1.95), 2.75, "*", width = 0.005, segmentParams = list(col = "#00468BFF"), textParams = list(col = "#00468BFF")) +
showSignificance( 2.3, c(min(d.sum.pptq$m) -0.01*min(d.sum.pptq$m),max(d.sum.pptq$m) + 0.001*max(d.sum.pptq$m)) , "*", width = -0.05)
#CDT Q
pcdtq <- plot_prepost_an(dm.cdt.quad.noout, "Cold detection threshold", "Log", df_an = an_mpp, v = "CDT_Quadriceps", x = F) + showSignificance( c(1.05,2.1), 0.37, "*", width = -0.005, segmentParams = list(col = "#ED0000FF"), textParams = list(col = "#ED0000FF")) +
showSignificance( 2.3, c(0.25, 0.3) , "*", width = -0.05)
# ADT
padt <- (plot_prepost_an(dm.adt, "Auditory detection threshold", "mA", df_an = an_mpp, v = "ADT_No site"))
p_mpp <- ppptq/pcdtq/padt +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(legend.position = "bottom",
plot.tag = element_text(family = 'Roboto', size = 11, face = "bold"))
ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/fig5_mpp.png", plot = p_mpp, dpi = 600, height = 10, width = 10)
p_mpp
4.4 Tables samples & models results
Results from models are summarized, along with summary statistics in a table. Data are shown using natural units. Two other tables are made to display results of log transformed CDT and PPT and of standardized values for each test_modality:site.
4.4.1 Original units
dt_mods <- data_long |>
select(ID, order, timing, condition, site, test_mod, test_value) |>
group_by(ID, order, site, test_mod) |>
pivot_wider(names_from = c(timing, condition), values_from = test_value) |>
rowwise() |>
mutate(baseline = pre_exp - pre_cont,
abc_exp = post_exp - pre_exp,
abc_cont = post_cont - pre_cont,
eih = (post_exp - pre_exp) - (post_cont - pre_cont)) |>
group_by(ID, site, test_mod) |>
pivot_longer(
cols = c(pre_exp, post_exp, abc_exp, pre_cont, post_cont, abc_cont, baseline, eih),
names_to = "cond",
values_to = "values"
) |>
group_by(ID, test_mod, site, cond) |>
summarise(ind_m = mean(values)) |>
group_by(test_mod, site, cond) |>
summarise(m = specify_decimal(mean(ind_m, na.rm = T), 2),
sd = specify_decimal(sd(ind_m, na.rm = T), 2)) |>
mutate(msd = glue("{m} ± {sd}")) |>
select(-m, -sd) |>
pivot_wider(names_from = cond, values_from = msd) |>
mutate(site = if_else(test_mod == "adt", "no site", site)) |>
mutate(test_mod = toupper(test_mod),
test_mod = factor(test_mod, levels = c("PPT", "PP", "HPT", "CDT", "ADT"))) |>
mutate(site = factor(site, levels = c("forearm", "quad", "no site"), labels =c("Forearm", "Quadriceps", "No site"))) |> arrange(site, test_mod)
p_int_mcs <- an_mcs |>
mutate(test_mod = factor(test_mod, levels = levels(dt_mods$test_mod))) |>
arrange(test_mod) |>
filter(term == "condition:site") |>
mutate(across(c(DenDF, statistic), round, 2)) |>
mutate(stat = glue("{statistic}~{NumDF},{DenDF}~")) |>
mutate(p.value = pvalue_format_table(p.value)) |>
select(test_mod, stat, p.value)
p_int_mpp <- an_mpp |>
mutate(test_mod = factor(test_mod, levels = levels(dt_mods$test_mod))) |>
arrange(test_mod) |>
filter(term == "timing:condition") |>
mutate(across(c(DenDF, statistic), round, 2)) |>
mutate(stat = glue("^#^{statistic}~{NumDF},{DenDF}~")) |>
mutate(p.value = pvalue_format_table(p.value)) |>
mutate(p.value = glue("^#^{p.value}")) |>
filter(test_mod == "ADT") |>
mutate(test_mod = factor(test_mod), site = factor(site)) |>
select(test_mod,site, stat, p.value)
p_int_mfull <- df_an |>
filter(effects == "test_mod:site") |>
rename(statistic = "F value") |>
rename(p.value = "Pr(>F)") |>
mutate(across(c(DenDF, statistic), round, 2)) |>
mutate(stat = glue("{statistic}~{NumDF},{DenDF}~")) |>
mutate(p.value = pvalue_format_table(p.value))
dt_mods <- full_join(dt_mods, p_int_mcs) |>
rename(fpp = stat) |>
rename(ppp = p.value)
dt_mods <- full_join(dt_mods, p_int_mpp) |>
mutate(
fpp = if_else(test_mod == "ADT",stat, fpp ),
ppp = if_else(test_mod == "ADT",p.value, ppp )
) |>
select(-stat, -p.value)
dt_mods$ffull <- c(rep(p_int_mfull$stat, 8), NA)
dt_mods$pfull <- c(rep(p_int_mfull$p.value, 8), NA)
t3 <- dt_mods |>
select(test_mod,site, pre_cont, post_cont, abc_cont, pre_exp, post_exp, abc_exp, fpp, ppp, baseline, eih, ffull, pfull) |>
arrange(test_mod) |>
flextable() |>
set_header_labels(
test_mod = "Test modalities",site = "Site", pre_cont = "Pre", post_cont = "Post" , abc_cont = "Diff~post-pre~", pre_exp = "Pre", post_exp = "Post" , abc_exp = "Diff~post-pre~", fpp = "F~df1,df2~", ppp = "p", baseline = "Baseline~diff~", eih = "EIH", ffull = "F~df1,df2~", pfull = "p"
) |>
add_header_row(values = c(
"Test modalities","Site", "Control", "Exercise", "Site x Condition interaction", "Baseline~diff~", "EIH", "Test modalities x Site interaction"
), colwidths = c(1,1,3,3,2,1,1,2)) |>
colformat_md(part = "all", j = 1:14) |>
merge_v(part = "header") |>
bold(part = "header") |>
align(part = "header", align = "center") |>
merge_v(j = c(1,2, 9,10,13,14), part = "body") |>
valign(j = c(1,2, 9,10,13,14), valign = "top") |>
hline(i = c(2,4,6), part = "body", border = dash_border) |>
hline(i = 8, part = "body") |>
vline(j = c(2, 10), border = dash_border, part = "body") |>
autofit() |>
set_table_properties(layout = "autofit") |>
add_footer_lines(values = c(
"Data from 38 participants.\n # F statistics and p value reported from timing x condition interaction. Units: CDT, |°C|; HPT, °C; PP, VAS 100 - 0; PPT, kPa.\n Abbreviations: CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Mechanical pinprick rating; PPT, Pressure pain threshold; Diff post - pre, Difference of post and pre values in each condition.") )
t3Test modalities | Site | Control | Exercise | Site x Condition interaction | Baselinediff | EIH | Test modalities x Site interaction | ||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Pre | Post | Diffpost-pre | Pre | Post | Diffpost-pre | Fdf1,df2 | p | Fdf1,df2 | p | ||||
PPT | Forearm | 441.25 ± 146.02 | 438.97 ± 138.12 | -2.28 ± 70.42 | 433.03 ± 131.34 | 440.49 ± 158.10 | 7.46 ± 101.14 | 9.921,414.09 | 0.00175 | -8.23 ± 107.34 | 9.75 ± 112.66 | 5.393,866.08 | 0.00111 |
Quadriceps | 680.39 ± 174.72 | 650.86 ± 182.37 | -29.54 ± 91.84 | 663.92 ± 183.60 | 709.87 ± 227.51 | 45.95 ± 101.21 | -16.47 ± 118.37 | 75.48 ± 130.82 | |||||
PP | Forearm | 63.85 ± 15.70 | 62.31 ± 17.63 | -1.55 ± 10.21 | 64.44 ± 16.22 | 62.52 ± 18.65 | -1.92 ± 11.09 | 01,414.09 | 0.949 | 0.58 ± 9.45 | -0.37 ± 12.82 | ||
Quadriceps | 70.38 ± 15.88 | 69.19 ± 15.27 | -1.19 ± 7.29 | 71.80 ± 14.40 | 69.82 ± 17.43 | -1.98 ± 10.70 | 1.42 ± 10.54 | -0.79 ± 12.18 | |||||
HPT | Forearm | 43.24 ± 2.45 | 43.00 ± 2.42 | -0.23 ± 1.72 | 42.55 ± 2.71 | 42.85 ± 2.53 | 0.30 ± 1.67 | 0.541,416.28 | 0.465 | -0.68 ± 1.94 | 0.53 ± 2.65 | ||
Quadriceps | 45.82 ± 2.86 | 45.58 ± 2.73 | -0.23 ± 1.72 | 45.41 ± 3.21 | 45.23 ± 2.91 | -0.17 ± 1.59 | -0.41 ± 1.58 | 0.06 ± 2.37 | |||||
CDT | Forearm | 1.43 ± 0.59 | 1.47 ± 0.63 | 0.04 ± 0.36 | 1.46 ± 0.72 | 1.52 ± 0.67 | 0.06 ± 0.58 | 7.61,414.55 | 0.00608 | 0.02 ± 0.44 | 0.02 ± 0.65 | ||
Quadriceps | 1.99 ± 0.92 | 2.10 ± 1.04 | 0.11 ± 0.67 | 2.19 ± 1.25 | 1.99 ± 1.13 | -0.21 ± 0.72 | 0.20 ± 0.82 | -0.32 ± 0.89 | |||||
ADT | No site | 0.35 ± 0.02 | 0.35 ± 0.01 | 0.00 ± 0.02 | 0.35 ± 0.01 | 0.35 ± 0.01 | 0.00 ± 0.01 | #1.111,412.44 | #0.294 | 0.00 ± 0.01 | 0.00 ± 0.02 | ||
Data from 38 participants. | |||||||||||||
#print(t2, preview = "docx")
save_as_docx(t3, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t2.docx")4.4.2 Log transformed CDT and PPT
PPT and CDT values were log transformed, values are reported in the table below.
t_log <- data_long |>
select(ID, order, timing, condition, site, test_mod, test_value) |>
ungroup() |>
rowwise() |>
mutate(test_value = if_else(test_mod %in% c("ppt", "cdt"), logt(test_value), test_value)) |>
group_by(ID, order, site, test_mod) |>
pivot_wider(names_from = c(timing, condition), values_from = test_value) |>
rowwise() |>
mutate(baseline = pre_exp - pre_cont,
abc_exp = post_exp - pre_exp,
abc_cont = post_cont - pre_cont,
eih = (post_exp - pre_exp) - (post_cont - pre_cont)) |>
group_by(ID, site, test_mod) |>
pivot_longer(
cols = c(pre_exp, post_exp, abc_exp, pre_cont, post_cont, abc_cont, baseline, eih),
names_to = "cond",
values_to = "values"
) |>
group_by(ID, test_mod, site, cond) |>
summarise(ind_m = mean(values)) |>
group_by(test_mod, site, cond) |>
summarise(m = specify_decimal(mean(ind_m, na.rm = T), 2),
sd = specify_decimal(sd(ind_m, na.rm = T), 2)) |>
mutate(msd = glue("{m} ± {sd}")) |>
select(-m, -sd) |>
pivot_wider(names_from = cond, values_from = msd) |>
mutate(site = if_else(test_mod == "adt", "no site", site)) |>
mutate(test_mod = toupper(test_mod),
test_mod = factor(test_mod, levels = c("PPT", "PP", "HPT", "CDT", "ADT"))) |>
mutate(site = factor(site, levels = c("forearm", "quad", "no site"), labels =c("Forearm", "Quadriceps", "No site"))) |> arrange(site, test_mod) |>
filter(test_mod %in% c("CDT", "PPT")) |>
select(test_mod,site, pre_cont, post_cont, pre_exp, post_exp) |>
arrange(test_mod) |>
flextable() |>
set_header_labels(
test_mod = "Test modalities",site = "Site", pre_cont = "Pre", post_cont = "Post" , pre_exp = "Pre", post_exp = "Post" ) |>
add_header_row(values = c(
"Test modalities","Site", "Control", "Exercise"
), colwidths = c(1,1,2,2)) |>
colformat_md(part = "all", j = 1:6) |>
merge_v(part = "header") |>
bold(part = "header") |>
align(part = "header", align = "center") |>
merge_v(j = c(1,2), part = "body") |>
valign(j = c(1,2), valign = "top") |>
hline(i = c(2), part = "body", border = dash_border) |>
vline(j = c(2), border = dash_border, part = "body") |>
autofit() |>
set_table_properties(layout = "autofit") |>
add_footer_lines(values = c(
"Data from 38 participants.\n Units: Log.\n Abbreviations: PPT, Pressure pain threshold; CDT, Cold detection threshold.") )
t_logTest modalities | Site | Control | Exercise | ||
|---|---|---|---|---|---|
Pre | Post | Pre | Post | ||
PPT | Forearm | 2.62 ± 0.15 | 2.62 ± 0.14 | 2.62 ± 0.12 | 2.61 ± 0.16 |
Quadriceps | 2.81 ± 0.12 | 2.79 ± 0.13 | 2.80 ± 0.12 | 2.83 ± 0.14 | |
CDT | Forearm | 0.15 ± 0.15 | 0.16 ± 0.15 | 0.15 ± 0.16 | 0.17 ± 0.16 |
Quadriceps | 0.27 ± 0.18 | 0.30 ± 0.18 | 0.31 ± 0.18 | 0.26 ± 0.21 | |
Data from 38 participants. | |||||
save_as_docx(t_log, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t2_log.docx")
#print(t2_log, preview = "docx")4.4.3 Standardized units
Models on EIH and absolute change used standardized values to take into account test modalities and site variability differences. PPT and CDT were log transformed prior to this procedure.
Values were standardized in each category of test modality and site using the following formula:
\[x_{standardized} = \frac{x_{obs} - \overline{x}}{SD}\] Leading each site:test modality level to have \(\overline{x} = 0\) and \(SD = 1\) as can be seen with the code below:
std_dat <- data_long |>
select(ID, order, timing, condition, site, test_mod, test_value) |>
ungroup() |>
rowwise() |>
mutate(test_value = if_else(test_mod %in% c("ppt", "cdt"), logt(test_value), test_value)) |>
group_by(test_mod, site) |>
mutate(test_value = scale(test_value)) # scale function standardized
std_dat |>
group_by(test_mod, site) |>
summarise(
m = mean(test_value, na.rm = T),
sd = sd(test_value, na.rm = T)
)# A tibble: 9 × 4
# Groups: test_mod [5]
test_mod site m sd
<fct> <fct> <dbl> <dbl>
1 ppt forearm -1.96e-14 1
2 ppt quad -5.12e-15 1
3 pp forearm 4.50e-15 1.00
4 pp quad 2.00e-15 1
5 hpt forearm 1.16e-14 1.00
6 hpt quad -4.26e-16 1.00
7 cdt forearm 1.32e-16 1.00
8 cdt quad 5.18e-16 1
9 adt no site -5.93e-14 1.00A third table of summary statistics is computed the pre, post, absolute change values in each condition and timing levels after standardization.
dt_mods_std <- std_dat |>
group_by(test_mod, timing, condition, site) |>
ungroup() |>
group_by(ID, order, site, test_mod) |>
pivot_wider(names_from = c(timing, condition), values_from = test_value) |>
rowwise() |>
mutate(baseline = pre_exp - pre_cont,
abc_exp = post_exp - pre_exp,
abc_cont = post_cont - pre_cont,
eih = (post_exp - pre_exp) - (post_cont - pre_cont)) |>
group_by(ID, site, test_mod) |>
pivot_longer(
cols = c(pre_exp, post_exp, abc_exp, pre_cont, post_cont, abc_cont, baseline, eih),
names_to = "cond",
values_to = "values"
) |>
group_by(ID, test_mod, site, cond) |>
summarise(ind_m = mean(values)) |>
group_by(test_mod, site, cond) |>
summarise(m = specify_decimal(mean(ind_m, na.rm = T), 2),
sd = specify_decimal(sd(ind_m, na.rm = T), 2)) |>
mutate(msd = glue("{m} ± {sd}")) |>
select(-m, -sd) |>
pivot_wider(names_from = cond, values_from = msd) |>
mutate(site = if_else(test_mod == "adt", "no site", site)) |>
mutate(test_mod = toupper(test_mod),
test_mod = factor(test_mod, levels = c("PPT", "PP", "HPT", "CDT", "ADT"))) |>
mutate(site = factor(site, levels = c("forearm", "quad", "no site"), labels =c("Forearm", "Quadriceps", "No site"))) |> arrange(site, test_mod)
p_int_mcs <- an_mcs |>
mutate(test_mod = factor(test_mod, levels = levels(dt_mods$test_mod))) |>
arrange(test_mod) |>
filter(term == "condition:site") |>
mutate(across(c(DenDF, statistic), round, 2)) |>
mutate(stat = glue("{statistic}~{NumDF},{DenDF}~")) |>
mutate(p.value = pvalue_format_table(p.value)) |>
select(test_mod, stat, p.value)
p_int_mpp <- an_mpp |>
mutate(test_mod = factor(test_mod, levels = levels(dt_mods$test_mod))) |>
arrange(test_mod) |>
filter(term == "timing:condition") |>
mutate(across(c(DenDF, statistic), round, 2)) |>
mutate(stat = glue("^#^{statistic}~{NumDF},{DenDF}~")) |>
mutate(p.value = pvalue_format_table(p.value)) |>
mutate(p.value = glue("^#^{p.value}")) |>
filter(test_mod == "ADT") |>
mutate(test_mod = factor(test_mod), site = factor(site)) |>
select(test_mod,site, stat, p.value)
p_int_mfull <- df_an |>
filter(effects == "test_mod:site") |>
rename(statistic = "F value") |>
rename(p.value = "Pr(>F)") |>
mutate(across(c(DenDF, statistic), round, 2)) |>
mutate(stat = glue("{statistic}~{NumDF},{DenDF}~")) |>
mutate(p.value = pvalue_format_table(p.value))
dt_mods_std <- full_join(dt_mods_std , p_int_mcs) |>
rename(fpp = stat) |>
rename(ppp = p.value)
dt_mods_std <- full_join(dt_mods_std , p_int_mpp) |>
mutate(
fpp = if_else(test_mod == "ADT",stat, fpp ),
ppp = if_else(test_mod == "ADT",p.value, ppp )
) |>
select(-stat, -p.value)
dt_mods_std$ffull <- c(rep(p_int_mfull$stat, 8), NA)
dt_mods_std$pfull <- c(rep(p_int_mfull$p.value, 8), NA)
t3_std <- dt_mods_std |>
select(test_mod,site, pre_cont, post_cont, abc_cont, pre_exp, post_exp, abc_exp, fpp, ppp, baseline, eih, ffull, pfull) |>
arrange(test_mod) |>
flextable() |>
set_header_labels(
test_mod = "Test modalities",site = "Site", pre_cont = "Pre", post_cont = "Post" , abc_cont = "Diff~post-pre~", pre_exp = "Pre", post_exp = "Post" , abc_exp = "Diff~post-pre~", fpp = "F~df1,df2~", ppp = "p", baseline = "Baseline~diff~", eih = "EIH", ffull = "F~df1,df2~", pfull = "p"
) |>
add_header_row(values = c(
"Test modalities","Site", "Control", "Exercise", "Site x Condition interaction", "Baseline~diff~", "EIH", "Test modalities x Site interaction"
), colwidths = c(1,1,3,3,2,1,1,2)) |>
colformat_md(part = "all", j = 1:14) |>
merge_v(part = "header") |>
bold(part = "header") |>
align(part = "header", align = "center") |>
merge_v(j = c(1,2, 9,10,13,14), part = "body") |>
valign(j = c(1,2, 9,10,13,14), valign = "top") |>
hline(i = c(2,4,6), part = "body", border = dash_border) |>
hline(i = 8, part = "body") |>
vline(j = c(2, 10), border = dash_border, part = "body") |>
autofit() |>
set_table_properties(layout = "autofit") |>
add_footer_lines(values = c(
"Data from 38 participants.\n # F statistics and p value reported from timing x condition interaction. \n Units: standardized; PPT and CDT were log transformed before standardization.\n Abbreviations: CDT, Cold detection threshold; HPT, Heat pain threshold; PP, Mechanical pinprick rating; PPT, Pressure pain threshold; Diff post - pre, Difference of post and pre values in each condition.") )
t3_std Test modalities | Site | Control | Exercise | Site x Condition interaction | Baselinediff | EIH | Test modalities x Site interaction | ||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
Pre | Post | Diffpost-pre | Pre | Post | Diffpost-pre | Fdf1,df2 | p | Fdf1,df2 | p | ||||
PPT | Forearm | 0.01 ± 0.99 | 0.02 ± 0.91 | 0.01 ± 0.51 | -0.01 ± 0.82 | -0.02 ± 1.03 | -0.01 ± 0.71 | 9.921,414.09 | 0.00175 | -0.01 ± 0.69 | -0.02 ± 0.81 | 5.393,866.08 | 0.00111 |
Quadriceps | 0.03 ± 0.93 | -0.12 ± 0.94 | -0.16 ± 0.51 | -0.04 ± 0.91 | 0.14 ± 1.02 | 0.18 ± 0.46 | -0.08 ± 0.63 | 0.34 ± 0.68 | |||||
PP | Forearm | 0.03 ± 0.87 | -0.05 ± 0.98 | -0.09 ± 0.57 | 0.06 ± 0.90 | -0.04 ± 1.04 | -0.11 ± 0.62 | 01,414.09 | 0.949 | 0.03 ± 0.53 | -0.02 ± 0.71 | ||
Quadriceps | 0.01 ± 0.93 | -0.06 ± 0.89 | -0.07 ± 0.43 | 0.09 ± 0.84 | -0.03 ± 1.02 | -0.12 ± 0.62 | 0.08 ± 0.61 | -0.05 ± 0.71 | |||||
HPT | Forearm | 0.12 ± 0.91 | 0.03 ± 0.90 | -0.09 ± 0.64 | -0.13 ± 1.01 | -0.02 ± 0.94 | 0.11 ± 0.62 | 0.541,416.28 | 0.465 | -0.26 ± 0.72 | 0.20 ± 0.99 | ||
Quadriceps | 0.10 ± 0.94 | 0.02 ± 0.90 | -0.08 ± 0.57 | -0.03 ± 1.06 | -0.09 ± 0.96 | -0.06 ± 0.52 | -0.14 ± 0.52 | 0.02 ± 0.78 | |||||
CDT | Forearm | -0.05 ± 0.89 | 0.01 ± 0.89 | 0.06 ± 0.55 | -0.03 ± 0.95 | 0.06 ± 0.94 | 0.09 ± 0.69 | 7.61,414.55 | 0.00608 | 0.02 ± 0.64 | 0.04 ± 0.87 | ||
Quadriceps | -0.05 ± 0.86 | 0.06 ± 0.89 | 0.11 ± 0.50 | 0.13 ± 0.87 | -0.14 ± 1.03 | -0.27 ± 0.68 | 0.18 ± 0.64 | -0.37 ± 0.77 | |||||
ADT | No site | 0.02 ± 1.21 | 0.09 ± 1.00 | 0.09 ± 1.57 | -0.01 ± 0.78 | -0.10 ± 0.71 | -0.08 ± 1.05 | #1.111,412.44 | #0.294 | -0.03 ± 0.75 | -0.16 ± 1.38 | ||
Data from 38 participants. | |||||||||||||
save_as_docx(t3_std, path = "/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/tables_wp1/t2_std.docx")5 Variables involved in EIH
An exploratory analysis is conducted to observe potential EIH contributors. As EIH was only observed at the quadriceps using PPT we consider the mean individual EIH (of three values) and relate it, using spearman rho correlation coefficient and scatterplots, to questionnaires scores, relative peak power output, rating of perceived exertion, changes in blood pressure, and individual data. To preserve statistical power for a purely exploratory analysis, no p value correction was applied.
# selecting data of personal characteristics
d_char <- data |>
select(ID, W_peak_rel, STAI_2, PCS, FPQ, IPAQ_total_moderate, IPAQ_total_vigourous, IPAQ_total_pa, TSK, bmi) |>
group_by(ID) |>
slice(1) |>
ungroup() |>
mutate(across(-ID, scale)) # scaling personal characteristics data
# selecting not controlled cardiovascular data and computing delta
d_cv <- d_cov_long |>
select(ID,condition, timing, br, sbp, dbp) |>
pivot_wider(names_from = c(condition, timing), values_from = c(br, sbp, dbp)) |>
mutate(
delta_br = (br_exp_post - br_exp_pre) - (br_cont_post - br_cont_pre),
delta_bp_diast = (dbp_exp_post - dbp_exp_pre) - (dbp_cont_post - dbp_cont_pre),
delta_bp_sys = (sbp_exp_post - sbp_exp_pre) - (sbp_cont_post - sbp_cont_pre),
) |>
select(ID, delta_br, delta_bp_diast, delta_bp_sys) |>
ungroup() |>
mutate(across(-ID, scale)) # scaling cardiovascular data
# select rpe and stai1
d_session <- data |>
select(ID, condition, RPE) |>
mutate(RPE = as.numeric(RPE)) |>
pivot_wider(names_from = condition, values_from = RPE) |>
mutate(delta_rpe = exp - cont) |>
select(ID, delta_rpe) |>
ungroup() |>
mutate(across(-ID, scale))# scaling rpe data
# a function is created to filter eih by site and test modality
# this data set is then merged with data sets from other variables
d_corr <- function(tm, s) {
d <- dm |>
filter(test_mod == tm, site == s) |>
group_by(ID) |>
summarise(
eih = mean(eih),
baseline = mean(test_value_pre),
delta_temp = mean(delta_temp)
) |>
select(ID, eih, baseline, delta_temp)
d1 <- full_join(d, d_char)
d2 <- full_join(d1, d_cv)
d_corr <- full_join(d2, d_session) |>
pivot_longer(cols = c(-eih, -ID), values_to = "val", names_to = "var2")
return(d_corr)
}
var_new_names_eih <- c( "W_peak_rel" = "Relative PPO" , "PCS" = "PCS total" , "FPQ" = "FPQ" , "IPAQ_total_moderate"= "IPAQ total mPA" , "IPAQ_total_vigourous" = "IPAQ total vPA" , "IPAQ_total_pa"= "IPAQ total PA", "W_peak" = "PPO", "TSK" = "TSK", "STAI_1" = "STAI-Y1", "bmi" = "BMI", "m_pre_quad_ppt" = "Baseline PPT", "m_pre_quad_cdt" = "Pre quad CDT", "m_pre_quad_hpt" = "Pre quad HPT", "m_pre_quad_pp" = "Pre quad PP", "m_pre_forearm_ppt" = "Pre forearm PPT", "m_pre_forearm_cdt" = "Pre forearm CDT", "m_pre_forearm_hpt" = "Pre forearm HPT", "m_pre_forearm_pp" = "Pre forearm PP","delta_rpe" = "Change RPE" , "delta_bp_sys" = "Change SBP", "delta_bp_diast" = "Change DBP" , "delta_temp" = "Change Tsk dominant quadriceps", "STAI_2" = "STAI Y2")cor.eihpptq <- d_corr("PPT", "Quadriceps") |>
pivot_wider(names_from = var2, values_from = val) |>
cor_test(
vars = "eih",
vars2 = c( "baseline","bmi","PCS", "FPQ", "TSK", "IPAQ_total_pa", "IPAQ_total_moderate", "IPAQ_total_vigourous", "W_peak_rel", "delta_bp_sys", "delta_bp_diast","delta_br", "delta_rpe", "delta_temp", "STAI_2"),
method = "spearman"
)
plot.corr.eihpptq <- full_join( d_corr("PPT", "Quadriceps"), cor.eihpptq) |>
filter(!var2 %in% c("bmi", "baseline", "delta_br")) |>
mutate(var2 = str_replace_all(var2, var_new_names_eih)) |>
filter(var2 != "Change Tsk dominant quadriceps") |>
mutate(var2 = factor(var2,
levels = c(
"Change SBP", "Change DBP", "Change RPE", "IPAQ total mPA", "IPAQ total vPA", "IPAQ total PA","Relative PPO", "FPQ", "PCS total", "STAI Y2", "TSK"
))) |>
ungroup() |>
rowwise() |>
mutate(
signif = pvalue_format_plot(p),
signif = factor(signif),
cor = round(cor, 2),
p = pvalue_format_plot(p),
corr = glue("rho = {cor}, {p}")) |>
ggplot(aes(y = eih, x = val)) +
geom_point(alpha = 0.7) +
geom_smooth(se = F, method = "lm", col = "#ADB6B6FF") +
geom_text(aes(label = corr, y = Inf, x = -Inf ), hjust = -0.1, vjust = +1.2, family = "Roboto", col = "black") +
labs(col = "P correlation", y = "Standardized EIH (PPT quadriceps)", x = "Standardized units") +
#scale_y_continuous(limits = c(-2,3)) +
facet_wrap(~var2, scales = "free") +
my_theme()
ggsave("/Users/vladimiraron/Library/CloudStorage/OneDrive-UCL/PhD/Projects/WP1/writing/figures_wp1/fig6_corr.png", plot = plot.corr.eihpptq, dpi = 1200, height = 8, width = 9)
plot.corr.eihpptq